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I. INTRODUCTION 


A. PROBLEM STATEMENT 


The dynamic response of a submersible vehicle operating at the extremes 
of its operational envelope is becoming increasingly important in order to en- 
hance vehicle operations. ‘Traditionally, dynamic stability of motion is stud- 
ied using eigenvalue analysis where the equations of motion are linearized 
around nominal straight line level flight paths (Arentzen & Mandel, 1960), 
(Clayton & Bishop, 1982), (Feldman, 1987). Under certain simplified as- 
sumptions, a simple criterion G, > 0 can be obtained where the stability 
index G, is function of the hydrodynamic coefficients in heave and pitch. 
Values for the stability index can be computed by, 


M.(Z, +m) 
aa aT (1) 


Gra 
This index is analogous to the familiar stability coefficient for horizontal 
plane maneuvering (Lewis, 1989) and can be thought of as a high speed ap- 
proximation where the effect of the metacentric restoring moment is minimal. 
If the value of G, is greater than zero, the vehicle is dynamically stable. As 
we point out in the next chapter though, this is only a sufficient, and rather 
conservative condition for stability. It is not a necessary condition in the 
sense that G, < 0 indicates instability at infinite forward speed. It is quite 


possible that at normal operating speeds the vehicle might be directionally 


stable. Furthermore, G, < 0 indicates a divergent loss of stability which is 
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quite uncommon in the vertical plane. Most modern submarines exhibit a 
flutter-like instability at high speed, which can not be analyzed using the 
above simplified index. Divergent motions may develop in combined six 
degrees of freedom (Papoulias et al 1993) and their occurrence can not be 


analyzed by a single stability index. 


B. OBJECTIVES AND OUTLINE 


In this work we examine the problem of stability of motion with controls 
fixed in the vertical plane, with particular emphasis on the mechanism of 
loss of stability of straight line motion. The surge equation is decoupled 
from heave/pitch through a perturbation series approach (Bender & Orszag, 
1978). It is shown that loss of stability occurs in the form of generic bifurca- 
tions to periodic solutions (Guckenheimer & Holmes, 1983). Taylor expan- 
sions and center manifold approximations are employed in order to isolate 
the main nonlinear terms that influence system response after the initial loss 
of stability (Hassard & Wan, 1978). Integral averaging is performed in order 
to combine the nonlinear terms into a design stability coefficient (Chow & 
Mallet—Paret, 1977). Special attention is paid to the study of the quadratic 
drag terms as they constitute some of the main nonlinear terms of the equa- 
tions of motion. This is a unique feature of the open loop problem. In similar 
studies of loss of stability under closed loop depth control (Bateman, 1993) 


it was found that the main damping mechanism is provided through the use 


of control susrafces. The difficulty associated with the nonsmoothness of the 
absolute value nonlinearities of the quadratic drag forces is dealt with by 
employing the concept of generalized gradient (Clarke, 1983). This has the 
advantage of keeping the linear terms constant, unlike the linear/cubic ap- 
proximation typically used in ship roll motion studies (Dalzell, 1978), where 
the linear damping coefficient is a function of the assumed amplitude of mo- 


tion. 


Vehicle modeling in this work follows standard notation (Gertler & Hagen, 
1976), (Smith et al 1978), and numerical results are presented for the DARPA 
SUBOFF model (Roddy, 1990) for which a set of hydrodynamic coefficients 
and geometric properties is available. Furthermore, the baseline vehicle is 
marginally stable with controls fixed under normal operating conditions and 
can serve as a prime example for the techniques described in this work. The 
model has been experimentally validated for angles of attack on the hull 
between +15 deg., while the constant coefficient approximation introduces 
very little error in time domain simulations (Tinker, 1978). Unless otherwise 
mentioned, all results in this work are presented in standard dimensionless 
form with respect to the vehicle length € = 4.26 m, and nominal forward 


speed U = 2.44 m/sec. All angular deflections are shown in degrees. 


II. PROBLEM FORMULATION 


A. EQUATIONS OF MOTION 


Assuming that vehicle motion is restricted in the vertical plane, the math- 
ematical model consists of the coupled nonlinear heave and pitch equations 
of motion. In a moving coordinate frame fixed at the vehicle’s geometrical 
center, Newton’s equations of motion for a port/starboard symmetric and 
neutrally buoyant vehicle are expressed in dimensionless form as follows, 


m(w — ug — zeq — wgq) = Zqt Zyw+ Zaq+ Zww 


nose 


—C'p Pr b(x)(w — 2q)|w — xq| dx + Z6 , (2) 
Ig + mzg(u + wq) — meg(w —uq) = Mjq+Myw+ M,q+ Mw 
+C'p a b(x)(w — 2q)|w — zq|z dz 


—zgpW cos 6 — zgpW sin 6 + Ms6 ; (3) 


where tgg = 2g — 2g, 2GB = 2G — Zp, and the rest of the symbols are based 
on standard notation and are explained in Table 1. Without loss of generality 
we can assume that zg = zg = 0, so that zgg = rg and zgg = 2g. The cross 
flow integral terms in these equations become very important for high angles 
of attack maneuvering, where they provide the primary motion damping. 
The drag coefficient, Cp, is assumed to be constant throughout the vehicle 


length for simplicity. This does not affect the qualitative properties of the 





results that follow. The vehicle pitch rate is, 
= q. (4) 


Dynamic coupling between surge and heave/pitch is present due to coordi- 
nate coupling as a result of the nonzero metacentric height. Therefore, pitch 


and heave motions must be studied together with surge, 


mu + mwq — mzgq + mzgq = 2 oad: + Xyu + Xwqgwg + 2 www? 


+Xyut’ = » ali X65 §° ’ (5) 


where we assume that both resistance and propulsive forces are proportional 


to the square of the speed or the propeller revolutions, respectively. 


In analyzing controls fixed stability of motion, the case 6 = 0 is examined 
first. The steady state solutions of the equations are determined by w = q = 
i = 9 = q = 0, where subscript 0 indicates variable value at steady state. 


Substituting these conditions in (2) we get, 
Zw Wo = C'p Ay Wo [evo | =) 5 (6) 


where, 
= nose d . 7 
Ay =f Wx) de (7) 
is the “waterplane” area. Since Z,, < 0, equation (7) admits only one solu- 
tion, namely wo = 0. Equation (3) then yields, 


tan 09 = GE : (8) 
2GB 


while (5) is used to determine the nominal forward speed, uo. 
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TABLE 1: NOMENCLATURE 


dummy independent variable 

steady state value of a 

expansion coefficients of zz in terms of z,, 22 
local beam of the hull 

quadratic drag coefficient 

regularization parameter 

stern plane deflection 

perturbation parameter, € = (mz2) 
criticality difference, € = u — u, 

vehicle mass moment of inertia 

nonlinear stability coefficient 

system eigenvalue 

vehicle mass 

pitch moment 

derivative of M with respect to a 

pitch rate 

polar coordinates of z,, 22 

transformation matrix of x to z 

pitch angle 

forward speed 

critical value of u 

heave velocity 

state variables vector, x = [0,w, q| 

surge force 

derivative of X with respect to a 

body fixed coordinates of vehicle center of buoyancy 
body fixed coordinates of vehicle center of gravity 
center of gravity/center of buoyancy separation, zg — rg 
vehicle metacentric height, z¢ — zg 

state variables vector in its normal form 
critical coordinates of z 

stable coordinate of z 

heave force 

derivative of Z with respect to a 

imaginary part of critical pair of eigenvalues 


> ns ee ee SS GS Se Me eS ee ee ee ee eg ee ee Ge eee ees eee 





B. REDUCTION OF ORDER 


The linearized surge, heave, and pitch equations of motion in the vicinity 


of the nominal point are, 


(m — X,)u + mzgq = 2Xyyu , (9) 
(m — Z,,)w — (mag + Z4;)=(Z,+m)q+ Zyw , (10) 
(I, — Mj)¢+mzgu —(mrzg + M,)w = Mw + (M, -— mze)q 
4+MQ6 , (11) 
where, 
Ms = topWi sin % — zggpW cos & , (12) 


is the hydrostatic restoring moment coefficient. The characteristic equation 


of (4), (9), (10), and (11) is obtained as, 


(—A,A + By )(A.r* + B.A? + C,A + Dz) + Agd*t + BzA* =0, (13) 
where, 
A; = m-—X,, 
B, = 741 GER ’ 


A, = (m— Z5)(1, — Mj) — (44 + mze)(Mu + mze), 

B, = —Zy(Iy — Mj) —(m— Zu)(M, — mag) — (Z, + m)(My + mzg) 
a MaiZ, + mre), 

See Mom — 25) + Zy(M, — mag) — M,(Z, +m) , 


Dy — Mo Zw ) 


Ag = (mzg)’(m— Zy), 


Bs = ~(mzg)Zy . 


It can be seen that the parameter (mzc) is responsible for surge and heave, 
pitch coupling. For zg = 0, equation (13) decouples into the surge eigenvalue 
A = B,/A, and the classical cubic characteristic equation for the vertical 
plane. It should be mentioned that the effect of the forward speed uw is 


embedded into the definition for the dimensionless vehicle weight W through, 


Ww — a . (14) 
If we introduce a smallness parameter, 
e=(mzc) , (15) 
we can rewrite (13) in the form, 
(A+ea)\*+(B+e8)\°+C0N+DA+4+E=0, (16) 


where in terms of previously defined coefficients, 


A = —Ajq@, 

B= —A,B,+B,A2 , 
C = -A,C,+ BiB , 
D = —-A,D,+ B,C? , 
E = Bids, 

a = m-Z, 

B= —Ly 








Following (Bender & Orszag, 1978) we expand the solutions of (16) in a 


regular perturbation series, 
A= Xo + AE + Ole”) ) (17) 


where Ao is an eigenvalue for ¢ = 0 (uncoupled surge or heave/pitch), , is 
the first order correction due to dynamic coupling, and O(e*) contains second 


and higher order terms in €. Substituting (17) into (16) we can get, 


NZ (aAo + B) 


i a 
2 AASE Bey sp a yop Ray, 


e+ O(e7). (18) 


It can be seen that the correction term is very small compared to the un- 
coupled root as evidenced by the \2 term. This is particularly true when 
Ao nears zero; i.e., close to a bifurcation point. Therefore, loss of stability 
can be studied by analyzing the heave/pitch equations decoupled from surge. 


The characteristic equation then becomes, 


fh ht SO Ge a ORD We 5 F = (0. (19) 


Plots of the system eigenvalues at nominal speed versus zg are shown in 
Figures 1 through 4. The surge eigenvalues is real and negative throughout 
the range of z¢, while the heave/pitch eigenvalues are real for small values 
of zg. The two larger real heave/pitch eigenvalues coalesce into a complex 
conjugate pair whose real part crosses zero for a certain value of zg. Within 
the accuracy of Figures 1 through 3, the eigenvalues A are identical to those 
computed by either the coupled or the uncoupled system, or the perturbation 


equations (18). There is a very small difference for the surge eigenvalue as 
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shown in Figure 4, but again the agreement between coupled and uncoupled 


systems is very good. 


C. DEGREE OF STABILITY 


The eigenvalues of the reduced order system (18) designate stability or 
instability of motion. A single measure of stability, the “degree of stability”, 
e,, can be defined as the maximum real part of all three eigenvalues of (18). 
This measures the slowest exponential convergence to the equilibrium when 
negative and the fastest exponential divergence from the equilibrium when 
positive. For all numerical calculations in this work, the degree of stability 
corresponds to the real part of a complex conjugate pair of eigenvalues. Typ- 
ical results are presented in Figures 5 through 7. Figure 5 shows the degree 
of stability versus zgg constant dimensionless speed up = 0.5 and different 
values of zgg. Similar results are shown in Figure 6 for constant zgg = 0.015 
and different values of uo. Finally, a three dimensional representation is 


depicted in Figure 7. 


D. CRITICAL SPEED 


The parameter value where the real part of the complex conjugate pair of 
eigenvalues shown in the previous figures crosses zero defines the point where 


linear stability is lost. This critical point can be computed by considering 
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Figure 1: System eigenvalues versus zG 
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Heave/pitch e-values 
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Figure 2: Complex conjugate heave/pitch eigenvalues versus z¢ 
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Figure 3: Real heave/pitch eigenvalues versus 2G 
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Figure 4: Surge eigenvalue versus zg 
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Figure 5: Degree of stability versus zgg for up = 0.5 and different values of z¢g 
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Figure 6: Degree of stability versus zgg for z¢g = 0.015 and different values of uo 
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Figure 7: Degree of stability versus zgg 
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equation (19). Routh’s criterion applied to this cubic yields A,D. = B.C, 


which can be solved for the dimensionless weight, 


B2C20 


—C 
Az D2, — B2C2, 


(20) 


where, 


Cr0 = Zw(M, = Mz) a 1 Oe + m) ; 
Co, = (m— Z,)(zGg cos 9% — zgp sin %) , 
Do, = Zw(eep sin % — zeg cos 4) . 


The value of the critical speed u, can then be evaluated from (20) and (14). 
Typical results are presented in Figures 8 and 9, nondimensionalized with 
respect to nominal vehicle speed 2.44 m/sec and length 4.26 m. Vertical plane 
motions are stable for forward speeds less than the critical speed. It can be 
seen that stability is increasing with increasing zg while zg = 0 is the most 
conservative condition for stability. Therefore, a vehicle which is stable when 
properly trimmed will remain stable for off—trim conditions. For comparison, 
we note that the simple stability coefficient G,, defined in equation (1), is 
monotonically decreasing and becomes more negative for decreasing zg, as 
shown in Figure 10. Thus it would have predicted unstable motions for the 
entire range of parameters shown in Figures 8 and 9. For completeness, the 
value of the steady state pitch angle, 6 (in degrees), is shown in Figure 11 
versus gg using zcg as the parameter. This pitch angle is computed from 


equation (8). 
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Figure 10: Stability coefficient G, versus zg¢g 
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Figure 11: Steady state pitch angle % versus zqg for z¢g from 0.005 to 0.025 with 
increments of 0.005 


II. BIFURCATION ANALYSIS 


A. INTRODUCTION 


In all cases of stability loss of the previous chapter, one pair of com- 
plex conjugate eigenvalues of the corresponding eigenvalue problem crosses 
transversally the imaginary axis. A situation like this in which a certain 
parameter is varied such that the real part of one pair of complex conjugate 
eigenvalues of the linearized system matrix crosses zero, results in the sys- 
en leaving its steady state in an oscillatory manner. This loss of stability 
is called Hopf bifurcation and generically occurs in either supercritical or 
subcritical form. In the supercritical case, stable limit cycles are generated 
after the nominal straight line motion loses its stability. The amplitudes of 
these limit cycles are continuously increasing as the parameter distance from 
its critical value is increased. For small values of this criticality distance the 
resulting limit cycle is of small amplitude and differs little from the initial 
nominal state. In the subcritical case, however, stable limit cycles are gener- 
ated before the nominal state loses its stability. Therefore, depending on the 
initial conditions it is possible to diverge away from the nominal straight line 
path and converge towards a limit cycle even before the nominal motion loses 
its stability. This means that in the subcritical Hopf bifurcation case the do- 
main of attraction of the nominal state is decreasing and in fact it shrinks 


to zero as the critical point is approached. Random external disturbances of 
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sufficient magnitude can throw the vehicle off to an oscillatory steady state 
even though the nominal state may still remain stable. After the nominal 
state becomes unstable, a discontinuous increase in the magnitude of motions 
is observed as there exist no simple stable nearby attractors for the vehicle 
to converge to. Distinction between these two qualitatively different types of 
bifurcation is, therefore, essential in design. The computational procedure 
requires higher order approximations in the equations of motion and is the 


subject of this chapter. 


B. THIRD ORDER EXPANSIONS 


The nonlinear heave/pitch equations of motion (2), (3), and (4) are writ- 


ten in the form, 


er (21) 
WwW = aw +ajoq + a13(2GB cos 8 + zgg sin 8) + d,(w,q) 

+c,(w,q) ; (22) 
G = aq, w + ao2q + Go3(2Gg cos 8 + zeg sin 8) + d,(w,q) 


+C (w, q) ’ (23) 


where 
D, = (m— Zy4)(1, — Mj) — (mae + 23)(meg + My) , 
a; Da = (ea eee 
a2D, = (1, — Mj)(m+ Z,) + (mae + Z5)(M, — mae) , 
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a3D, = —(mzg+ Z;)W , 
aD, = (m—-Z,)M, + (meg + My)Zy , 
a.9D, = (m—- Z4)(M, — mag) + (mag + My)(m+ Z,) , 
oe (in — 25) , 
do(w,4)De = (ly — Mj)lw + (nee + Zale 
ea,g)D, = (m— Zy)1,+ (mae + Ma )l, , 
q(w,q)D, = (1, — Mj)mzeq’ — (mag + Z;)mzgugq , 
oo(w,q)D, = —(m— Zy)mzgwq + (meg + My)mzeq’ , 


and J,,, J, are the cross flow integrals 


Ly 


Cp [4 2)(w — 2q)|w — 2q| dz , (24) 


Cp ~— bx)(w — xq)|\w — zq|zxdz . (25) 


I, 


The system of equations (21) through (23) is written in the compact form 
x = Ax + g(x), (26) 
where 

es [ey ine| (27) 
is the three state variables vector, and A is the linearized sytem matrix eval- 
uated at the nominal point x9. The term g(x) contains all nonlinear terms 
of the equations. Hopf bifurcation analysis can be performed by isolating 
the primary nonlinear terms in g(x). Keeping terms up to third order, we 

can write 
gio) = g(x) + g(x). (28) 
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Using equations (21) through (25), the various terms in (28) can be written 


as, 
qe = 0 } 
gs”) — Gh — M;)mzeq7 — (mzg + Zy Nzgwq + d\?)(w,q) , (29) 
gs) = —(m—Zs)mzgwq + (mag + My)mzg@ + d?)(w,q) , 
and 
(a a0r 
gf = d\?) (w, q) a 4 a13(2GB sin 0) — ZgB cos 9 0° ’ (30) 
gs) = d)(w,q) + ta23(rep sin > — 2Gp Cos O)0° 


Expansion in Taylor series of d,, d, requires expansion of the cross flow 


integrals [,,, [,, which require the Taylor series of 


f(E) = €l€| . (31) 


This expression can be converted into an analytic function using Dalzell’s 


approximation (Dalzell, 1978), 
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i) 
El€| a 16 °°¢ = 48 £. ’ (32) 


which is derived by a least squares fit of an odd series over some assumed 
range of €, namely —€,. < € < €. This approximation, which is shown in 
Figure 12, has been extensively used in ship roll motion studies and is very 
useful for its intended purpose. However, in the present problem it suffers 


from several major drawbacks: 
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e It introduces a linear term which depends on the assumed range of 


motion, and it renders the critical speed function of the vehicle motions. 


e The cubic term, which is ultimately responsible for the Hopf bifurcation 
analysis, is a function of the assumed range of vehicle motions which 


can not be known in advance. 


e As Figure 12 demonstrates, the slope of the actual curve at the origin 
is significanly different than the approximation, which would make the 


bifurcation results unreliable. 


Instead of Dalzell’s approximation, we employ the concept of generalized 
gradient (Clarke, 1983), which is used in the study of control systems in- 
volving discontinuous or non-smooth functions. In this way we approximate 
the gradient of a non-smooth function at a discontinuity by a map equal to 
the convex closure of the limiting gradients near the discontinuity. In our 


problem we write, 


f(£) = €oléo| + 2lgol(é — £0) + sign(£o)(€ — £0)” + FE) , (33) 


as the Taylor series epansion of f(€) near &. The sign function in (33) can 


be approximated by, 
sign()) = im tanh (®) ; (34) 


A graphical representation of the approximation (34) is shown in Figure 


13. The quantity y is a small regularization parameter and is used in the 
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Figure 12: Graphical representation of Dalzell’s approximation of ¢|€| versus &/€. 
Solid curve is the exact expression and dotted curve is the approximation (32) 
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next section for the proper normalization of the results. Using (34), we can 


approximate f(£) in the vicinity of  =0 by, 
ele] — 2 (35) 
=e 


Since 


frw—2q, (36) 


we can express the non-smooth cross flow integral terms by, 


Cp 


= rm —(Eow® — 3E,w’q + 3E,wq? — Esq’) , (37) 
Cp 
= rm —(E\w*® — 3E,w’q + 3E3wq’? — Exq’) , (38) 


where 


E; = = [aif ‘p(z) de , (39) 


are the moments of the vehicle “waterplane” area. 


C. COORDINATE TRANSFORMATIONS 


Using the previous second and third order Taylor series expansions, equa- 


tion (26) is written in the form, 
Se Ax + g(x) + (0). (40) 


If T is the matrix of eigenvectors of A evaluated at the critical point u = u,, 


the linear change of coordinates, 


see se (41) 
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Figure 13: Graphical representation of the sign function and its hyperbolic tangent 
approximation (34). Solid curve corresponds to y = 0.1 and dotted curve corre- 
sponds to y = 0.01 


30 


transforms system (40) into its normal coordinate form, 
ae ng (Tz) PT (Tz). (42) 


At the Hopf bifurcation point, matrix T~' AT takes the form, 


0 —Wp 0 
TT AT = Wy 0 0 ; 
0 QO p 


where wo is the imaginary part of the critical pair of eigenvalues, and the 
remaining eigenvalue p is negative. For values of u close to the bifurcation 


poit u., matrix T~'AT becomes, 


a’ € —(wo + w’e) 0 
TAT = | (wo +u’e) a’ € 0 : 
0 0 p+pe 


where € denotes the criticality difference 


E=U—4U,, (43) 
and 
a’ = derivative of the real part of the critical. eigenvalue 
with respect to €, 
w' = derivative of the imaginary part of the critical eigenvalue 
with respect to €, 
p = derivative of p with respect toe. 
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D. CENTER MANIFOLD EXPANSIONS 


Due to continuity, the eigevalue p+ p’e remains negative for small nonzero 
values of «. Therefore, the coordinate z3 corresponds to a negative eigenvalue 
and is asymptotically stable. Center manifold theory predicts that the rela- 
tionship between the critical coordinates z,, z2 and the stable coordinate 2, 


is at least of quadratic order. We can then write 23 as, 
Zy = 11 ze + 122122 + Ql22 25 : (44) 


where the coefficients, a;;, in the quadratic center manifold expansion (44) 


need to be determined. By differentiating equation (44) we obtain, 


zg = 20412121 + 12(2Z122 + 222) + 20222222 « (45) 
We substitute z; = —wo22 and Zz. = woz, from equation (42) into (45), and 
we obtain 

23 = Ql, 2Wp 2? + 2(a22 — 011 \Wo 2 22 0119W9 23 ‘ (46) 


The third equation of (42) is written as, 


éy = pzy + [T-*g)(Tz)| (47) 


(3,3) ’ 


where terms up to second order have been kept. If we denote the elements 


of T and T™ by, 


T=(m;), T= [nj] , (48) 
then 
dy 
T-'g)(Tz) =} d, | , 
d3 
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where 


dy = m2(l2521 + L262 22 + bo725) + M13(la52; + l3621 22 + bs72z5) , (49) 
dy = 12(lo521 + lo621 22 + €2723) + N23(bas zi + b36 21 22 + re , (50) 


d3 = naz (log Zz: + bog 21 22 + lo7 25 ) aa m3 (£55 2, oe b36 21 22 + l37 25 ) : (51) 


Expressions for the coefficients £; are given in the Fortran programs in the 


appendix. 
Equation (47) then becomes 
23 = p23 Fr da , (52) 
and substituting (44) and (51) into (52) we get, 


23 = (pay; + Ngoblo5 + N33 35 Ve + (pare + N32log + 133 l36 )Z1 22 


+( parr + 132 £27 + 133 £37) z5 (53) 
Comparing coefficients of (46) and (53) we get, 


Pay, + Wo. = Ngzlo5 + Naglgs , (54) 
—2w 911 — Pay. + 2wy 22 = Nagle + Na3lse ’ (55) 


132627 + N33637 . (56) 


| 


—W9 12 — PX&22 


Solution of the system of linear equations (54) through (56) yields the coef- 


ficients in the center manifold expansion (44). 
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E. AVERAGING 


Using the previous Taylor expansions and center manifold approxima- 


tions, we can write the reduced two-dimensional system that describes the 


center manifold flow of (42) in the form, 


Zz = Wez, — (wo +w'e)z. + F(z, 22) , 
Z. = (wo tw'e)z, + a’ eza + Fo(z1, 22) , 
where 
FO air) T1124 + 11221 29 + 11321 24 + 1142) 
+pirz, + Pioziz2 + Pis2y , 
F,(21,22) = 12, ss 122.2) 29 is 1932124 =F T2425 
+pa, zi + poaz1Z2 + Pra2s 
and 


Tis = Nilo; + r43 43; , 1 = eee ,=1,. 2. 


Pij = Nirl + nisls,, t=1,2, g=1,2,3, k=7+4. 


If we introduce polar coordinates in the form, 


2, =Kcosd, 23.7 sin oe 


(57) 


(58) 


(59) 


(60) 


(61) 


(62) 


(63) 


we can use (57) and (58) to produce an equation describing the rate of change 


of the radial coordinate R, 


R=a'cR + P(d)R° + Q(4)R? . 
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(64) 


This equation contains one variable, R, which is slowly varying in time, and 
another variable, ¢, which is a fast variable. Therefore, equation (64) can be 
averaged over one complete cycle in ¢ to produce an equation with constant 


coefficients and similar stability properties, 
R=a'eR+ KR°+ LDR’, (65) 


where 


1 2n 
ce = = | P(¢) do = 5(3r 11 +1ry3 +792 + 3ro4) , (66) 


L a d 
= =| A¢)db =0. (67) 
Therefore, the averaged equation (65) becomes 


R=a'eR+KR’*. (68) 


F. LIMIT CYCLE ANALYSIS 


Equation (68) admits two steady state solutions, one at R = 0 which 


corresponds to the trivial equilibrium solution at zero, and one at 


Ro = ae < (69) 


This equilibrium solution corresponds to a periodic solution or limit cycle in 
the cartesian coordinates z,, 2). For this limit cycle to exist, the quantity Ro 
must be a real number. In our case a’ is always positive, since the system 


loses its stability; i.e., the real part of the critical pair of eigenvalues changes 
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from negative to positive, for increasing u. Therefore, existence of these 


periodic solutions depends on the value of K. Specifically, 


e if K <0, periodic solutions exist for « > 0 or u > u,, and 


e if K > 0, periodic solutions exist for « < 0 or u < u,. 


The characteristic root of (68) in the vicinity of (69) is 
3 =—2a’e, (70) 


and we can see that 


e if periodic solutions exist for u > u, they are stable, and 
e if periodic solutions exist for u < u, they are unstable. 
The period of these periodic solutions can be estimated as follows. Equa- 
tions (57), (58), and (63) produce an equation in ¢ similar to (64), 
$= uy +w'e+ F()R? + G(A)R. (71) 
The averaged form of (71) is 
6 = wy +w'e+ MR’, (72) 


where 
2n 
Vo a F(¢)d¢ = = (sri + 723 — T12 — 3ry4) - (73) 


The limit cycle period can be computed by substituting (69) into (72), 


20 20 wi K —a'M 
ee pe alae 74 
wo twet+t MRR wu ( wo ik ‘| FO tA) 
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G. RESULTS AND DISCUSSION 


Typical results of the nonlinear stability coefficient K are shown in Fig- 
ures 14 and 15. Figure 14 presents a plot of K -y versus zg for zg = 0.015 
and for different values of Cp. It should be emphasized that the use of K -¥ 
is more meaningful than the use of K, since it properly accounts for the 
use of the regularization parameter y as seen from equations (35) and (68). 
Numerical evidence demonstrates that all curves K -y versus zg converge 
for y — 0. For practical purposes, values of y smaller than 0.001 produce 
identical results. The results of Figure 14 demonstrate the profound effect 
that the quadratic drag coefficient Cp has on stability of limit cycles. All 
Hopf bifurcations are supercritical (K < 0), and they become stronger su- 
percritical as Cp is increased. It is worth noting that results for Cp = 0 
produce subcritical behavior, K > 0, which is clearly incorrect. Thus, ne- 
glecting the effects of Cp would have produce entirely wrong results in the 
present problem. Figure 15 shows a plot of K -+y versus zg for Cp = 0.5 
and different values of the metacentric height zg. It can be seen that, the 


bifurcations become stronger supercritical as initial stability z¢ is increased. 


The bifurcation analysis results are verified by direct numerical simula- 
tions shown in Figures 16 through 18. Figure 16 shows the results of two 
numerical simulations for two values of nominal speed up in terms of the 
vehicle pitch angle 6 (in degrees) versus time (in seconds). The critical value 
of speed, u,, is about 0.495 as can be seen from Figure 8, while the other 


parameters in the simulations were Cp = 0.5, zg = 0.015, and zg = 0. It 
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Figure 14: K -y versus zg for z = 0.015 and different values of the drag coefficient 
Cp 
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x10-3 


Figure 15: K -+ versus zg for Cp = 0.5 and different values of the metacentric 
height z¢ 
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Figure 16: Time histories (0,t) for Cp = 0.5, zg = 0.015, zq = 0, and two different 
values of nominal speed 
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Figure 17: Ti 
nominal forward speeds 
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Figure 18: Limit cycle amplitudes from the results of Figure 17 
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can be seen that convergence to zero is ensured for up < u, and convergence 
to a limit cycle occurs for uw > u,. This indicates supercritical behavior as 
shown before. A selection of time histories is shown in Figure 17 for a range 
of forward speeds and the same parameters as in Figure 16. The same initial 
disturbance, 9 = 5 degrees, was introduced at t = 0 for all simulations. It 
can be seen that the amplitude of limit cycles increases as the distance of u 
from u, is increasing. The rate of convergence of solutions to their limit cy- 
cles is also increasing, while their period remains essentially constant. These 
results are summarized in Figure 18, where the amplitudes of the numeri- 
cally computed limit cycles are plotted versus ug. The behavior is clearly 


) 
| 
supercritical, which agrees with our findings of the bifurcation analysis. 
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IV. BIAS EFFECTS 


A. LOSS OF STABILITY 


Stability analysis of motions at a nonzero angle of attack can be per- 
formed by first introducing some bias into the steady state solution and its 
perturbations. This can be achieved by maintaining a nonzero dive plane 
angle at nominal. In this case the steady state solutions are gg = 0, and wp, 


6) are computed from, 
Ly Wo = Cp Eowo|wo| + Z56 = 0 ; (75) 
Muwo + Cp FE, wolwo| — W(zeg cos 9 + zg sin %) + Msd = 0. (76) 


The coefficients Ey, E, are computed from (39). In order to solve (75) we 


observe that when 6 > 0, then wo < 0, 


—Zy = 4) Fee Ee 
D#046 (77) 


—2C p Eo 


Wo = 


and when 6 < 0, then wo > 0, 


7, =a) ZoeaO en aes 
D+046 . (78) 


—2C'p Eo 


Wo = 
The angle 9%) can then be computed from equation (76). 


Linearization of the equations of motion in the neighborhood of the above 


equilibrium point produces the linear system, 


(m — Zy3)w — (mag + Z3)¢ = (Jy — 2Cp Ep|wo|)w 


4b 





+(Z, +m + 2CpF, |wol)¢ , (79) 
—(My + mag)w + Uy — Mj)q = (My + 2Cp & |wo|)w 

+(M, — mtg — mzguo — 2Cp E|wol)¢ 

+W (aes sin 9 — zep cos 9 )6 , (80) 


=¢, (81) 


where the variables w, @ are understood as small deviations from their equi- 
librium values. Numerical solution of the generalized eigenvalue problem (79) 
through (81) yields the critical speed values where the nominal equilibrium 


solution becomes unstable. 


B. ANALYSIS OF HOPF BIFURCATIONS 


It can be numerically verified that the above calculations for the new 
critical speed result in a loss of stability in the form of Hopf bifurcations, as 
for the 6 = 0 case. These Hopf bifurcations can be analyzed using the same 
general methodology that was developed in the previous chapter. The non- 
linear expansions are like equations (21) to (23) with the following changes: 


The substitutions 
Lw — Ly —2WpEoro , 
M, ae M, a 2C'p E, wo ? 
Za +m — Zag +m+2Cp ku» , 
M,-—m«g — M,—-—mzg —mzguo — 2Cpkrr» , 
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are assumed in the definition of coefficients a;;. Furthermore, equation (33) 
prduces 2nd order contributions due to wo # 0, as well as 3rd order. Using 


(33) we can compute the 2nd order expansions of di?) and d”) in (29) using, 


I?) = Cp(Eow? —2E,wq + Ex’) , (82) 
I” 


Cp(E,w* — 2E,wq + E3q’) (83) 


These equations are valid for wo > 0 or 6 < 0. For 6 > O the signs of £; 
must be switched. The third order expansions J‘%) and I\3) are the same 
as in equations (37) and (38). Using these additional terms, the nonlinear 
stability coefficient AK can be computed in the same way as in the previous 


chapter. 


C. RESULTS AND DISCUSSION 


Typical results for 6 # 0 are shown in Figures 19 through 24. Figure 
19 shows the equilibrium pitch angle 4 (in degrees) versus zg for different 
values of 6 from —5 to 5 degrees with increments of 1 degree. Solid curves 
correspond to positive 6 and dashed curves to negative. Figure 20 shows the 
degree of stability for the equilibrium points of Figure 19, while Figure 21 
presents the degree of stability in a three dimensional view. It can be seen 
that positive and negative values of 6 have almost identical stability charac- 
teristics. Furthermore, the degree of stability becomes more negative as the 
absolute value of 6 is increased, which means that we expect a wider domain 


of stability in this case. This is verified by the critical speed plots shown in 
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Figure 20: Degree of stability versus zg for z = 0.015, uw = 0.5, and different 
values of 6 
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Figure 21: Degree of stability versus zg and uo for z@ = 0.015 and two values of 6 
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Figure 22: Critical speed u. versus zg for zg = 0.015 and different values of 6 
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Figure 23: Critical speed u. versus zg and 6 for z@ = 0.015 
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Figure 24: Nonlinear stability coefficient K versus zg for Cp = 0.5, zg = 0.015, 
and different values of 6 
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Figures 22 and 23. The critical speed is minimum for 6 = 0 and it increases 
monotonically with increasing absolute value of 6. This stabilizing effect of 
asymmetry (bias) remains approximately true for the nonlinear analysis, as 
demonstrated by the results of Figure 24. It can be seen that the nonlinear 
stability coefficient AK becomes more negative as 6 is decreasing from zero. 
For increasing 6, kK becomes less negative but the difference from the 6 = 0 
calculations appears to be very small. Therefore, limit cycle stability is not 
significantly affected by the bias effects that are induced by small nonzero 


dive plane angles. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


This thesis presented a comprehensive nonlinear study of straight line stabil- 
ity of motion of submersibles in the dive plane under open loop conditions. 
A systematic perturbation analysis demonstrated that the effects of surge in 
heave/pitch are small and can be neglected. Primary loss of stability was 
shown to occur in the form of Hopf bifurcations to periodic solutions. The 
critical speed were instability occurs was computed in terms of metacentric 
height, longitudinal separation of the centers of buoyancy and gravity, and 
the dive plane angle. Analysis of the periodic solutions that resulted from the 
Hopf bifurcations was accomplished through Taylor expansions, up to third 
order, of the equations of motion. A consistent approximation, utilizing the 
generalized gradient, was used to study the non—analytic quadratic cross flow 
integral drag terms. The results indicated that loss of stability occurs always 
in the form of supercritical Hopf bifurcations with stable limit cycles. It was 
shown that this is mainly due to the stabilizing effect of the drag forces at 
high angles of attack. As a recommendation for further research, we suggest 


a nonlinear analysis of coupled motions in six degrees of freedom. 


APPENDIX 


The following is a list and description of the computer programs used in this 


thesis. The programs are written in FORTRAN or MATLAB. Complete 


printouts of the programs follow after the list. 


PERTURB.M 
MATLAB program for performing surge/heave/pitch perturbation anal- 


ysis. 


CRIT_0.M 


MATLAB program for calculating the critical speed for 6 = 0. 


CRIT _DELTA.M 


MATLAB program for calculating the critical speed for 6 + 0. 


HOPF_0.FOR 
FORTRAN program for calculating the nonlinear stability coefficient 


iioro — |. 


HOPF _DELTA.FOR 
FORTRAN program for calculating the nonlinear stability coefficient 
K for 6 £ 0. 


SIM.FOR 


FORTRAN program for simulating the equations of motion. 


4) 


i, Program perturb.m 

7, This file applies concepts of perturbation analysis theory, 
4 in order to prove that the 3 by 3 system describing motion 
74 in the vertical plane, decouples from the 4 by 4 one. 


rho = 1.94; 

g = 3272; 

E = 13.9792; 

ndi = 0.6*rho*L"2; 

nd2 = 0.5*rho*L’ 3; 

nd3 = 0.5*rho*L"4; 

nd4 = 0.5*rho*L"§; 

iy = 661.32/nd4; 

zqdot = -6.330-4; 

zwdot = -1.4629e-2; 

zq = 7.5450-3; 

zw = -1.3910-2; 

mqdot = -8.8e-4; 

mwdot = -&.61¢e-4; 

mq = -3.7020e-3; 

mw = 1.0324e-2; 

cd = 0.015; 

zb = C/L; 

m = 1656.2363/(g*nd2) ; 
xudot = -0.05*m; 

xb = O/L; 

uo = 4; 

W = 1856.2363./(nd1.*uo.~2); 
b = wW; 

Xg 0; 

zg = linspace(0.001,0.025,500) ; 
xgb = xg-xb; 

zgb = zg-zb; 

theta = atan(-xgb./zgb); 
alpha = m-zwdot; 

beta = -Zw; 

Ail = m-xudot; 

B1 = -2*cd; 
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> 
I) 


(m-zwdot )*(iy-mqdot ) -(m*xg+zqdot )*(mwdot+m*xg) ; 
B = -zw*(iy-mqdot)-(m-zwdot )*(mq-m*xg)-(zq+m)*(mwdot+m*xg)... 
-(m*xg+zqdot) *mw; 


for i = 1:length(zg) 
deltaz(i) = xgb*w*sin(theta(i))-zgb(i)*w*cos(theta(i)); 
C(i) = -deltaz(i)*(m-zwdot)+zw* (mq-m*xg) -(zq+m) *mw; 
D(i) = deltaz(i)*zw; 
A2(i) = (m*zg(i))°2*(m-zudot) ; 
B2(i) = -(m*zg(i))~2*zw; 


CA = -A1*A; 

CB = -A1*B+Bix*A; 

CC = -Ai*C(i)+B1*B; 

CD = -A1*D(i)+B1i*C(i); 
CE Bi*xD(i); 


p (CA CB CC CD CE]; 
r(i:4,i) = roots(p); 


lamda(i:4,i) = r(i:4,i)-(r(1:4,i).°3.*(alpha.+*r(1:4,i)+beta))... 


(Cm. *zg(i)).72./(4.*CA.*r(1:4,1).73+... 
3.*CB.*r(1:4,1).°2+2.*CC.*r(1:4,1)+CD); 


PA = -A1*A+A2(i); 

PB = -A1*B+B1¥*A+B2(i) ; 
PC = -Ai*C(i)+Bix*B; 

PD = -A1*D(i)+B1*C(i) ; 


PE = Bi*D(i); 
p4 = [PA PB PC PD PE]; /, e-values of 4 by 4 system. 
r4(1:4,i) = roots(p4); 

end 


o¢ 


/, Program crit_O.m 
7, Evaluation of critical speed for delta=0 


rho = 1.94; 

g = 32.2; 

If = 13.9792; 

ndi = 0.5*rho*L~2; 
nd2 = 0.5*rho*L"3; 
nd3 = 0.5*rho*L"4; 
nd4 = 0.5*rho*L5; 


iy = 561.32/nd4; 
zqdot = -6.330e-4; 
zwdot = -1.4529e-2; 
zq = 7.545e-3; 
zw = -1.391e-2; 
mqdot = -8.8e-4; 
mwdot = -5.610e-4; 


mq = -3.702e-3; 

mw = 1.0324e-2; 

cd = 0.015; 

zb = O/L; 

m = 1666.2363/(g*nd2) ; 
xudot = -0.05*m; 

xb = O/L; 


xg = linspace(-0.01,0.01,40); 
Zg linspace(0.005,0.025,40); 


Gv = 1 - mw.*(zqtm)./(zw.*(mq-m.*xg)) ; 
lambda = -2*cd/m-xudot; 7% lambda = -1.6439 


xgb 
zgb 


xg-xb; 
zgZ-Zb; 


for j=i:length(zg) 
for i=1:length(xg) 


theta(i,j) = atan(-xgb(i)/zgb(j)); 

aQ = (m-zwdot)*(iy-mqdot )-(mwdot+m*xg (i) )*(zqdot+m*xg(i)); 

bO (-zwdot*m-m*mw-zq*m) *xg(i)+(-m*mq+zwdot*mq-zqdot*mw... 
-zqemudot-m*mwdot-iy*zwtmqdot*zw) ; 
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cO = -m*zw*xg(i)+mq*zw-zq*mw-m*mw ; 

cl = (-m*xg(i)+zwdot*xg(i)+m*xb-zwdot*xb)*sin(theta(i,j))... 
+(-m*zb-zwdot*zg(j)+tzwdot*zb+m*zg(j))*cos(theta(i,j)); 

di = (zw*xg(i)-zw*xb)*sin(theta(i,j))+... 


(zw*zb-zw*zg(j))*cos(theta(i,j)); 
’, After applying AD=BC... 


wOi,j) = bO*c0/(a0*d1-b0*c1) ; 
UOC ij) = (1556.2363/(ndi*w(i,j)))7~.5; 


end 
end 
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4 File crit_delta.m 
4, Evaluation of critical speed for nonzero stern plane angle. 


rho = 1.94; 

gz = 32.2; 

L = 13.9792; 

ndi = 0.6*rho*L~2; 
nd2 = 0.6*rho*L73; 
nd3 = 0.5*rho*L~4; 
nd4 = 0.5*rho*L5; 


iy = 661.32/nd4; 
zqdot = -6.330e-4; 
zwdot = -1.4529e-2; 
Zq = 7.6466-3; 
Zw = -1.3910e-2; 
mqdot = -8.8e6-4; 
mydot = -6.61e-4; 


mq = -3.7026-3; 

mw = 1.03246-2; 

m = 1656.2363/(g*nd2) ; 

xudot = -0.05*m; 

xb = O/L; 

zb = O/L; 

ZE = 0.016; 

cdz wee 0.6; 4 if cdz=0, division by zero provides NaN... 
zd = -65.6036e-03; 

md = -2.4096e-03; 

EO = 0.1036509; /, EO == Integral of b(x)*dx 

E1 = -2.36299826-03 ; 4 E1 == Integral of x*b(x)*dx 

E2 = 7.3672147e-03; , E2 == Integral of x°2*b(x)*dx 

xg = linspace(-0.01,0.01,80); 4 NON-DIMENSIONAL 

uO = linspace(2,5,40); 

xgb = xg-xb; 

zgb = zg-zb; 

wi = 1556.2363./(nd1.*u0.~2); 

delta = input(’Enter the value of delta: delta = linspace( , , )’); 


delta = delta*pi/180; 
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for k=1:length(delta) 
for j=1:length(u0) 
for i=1:length(xg) 


if delta(k) > 0 
wO = -(zwtsqrt (zw 2-4*cdz*EO*zd*delta(k)))/(2*cdz*E0) ; 
coeffi = (mw*wO-cdz*w0” 2*Ei+md*delta(k))/wi(j); 
poly1 ((xgb(i)“2+zgb°2) (-2*coeffi*zgb) (coeffi ~2-xgb(i)~2)]; 
thetaO = asin(roots(poly1)); 
checki = xgb(i)*cos(theta0(1))+zgb*sin(theta0(1)); 
check2 = xgb(i)*cos(theta0(2))+zgb*sin(theta0(2)); 
if abs(check2-coeff1) < abs(checki-coeff1) 
thetaO = theta0(2) ; 
else 
thetaO = theta0(1); 
end 


B=[ (m-zwdot) -(m*xg(i)+zqdot) 0 
-(mwdot+m*xg(i)) iy-mqdot 0 
0 0 1]; 


A= ((2zw+2*cdz*E0*w0) (zqt+m) -2*cdz*E1*w0 0 
mw-2*cdz*E1*wO (mq-m*xg(i))-m*zg*wO+2*cdz*E2*wO ... 
(xgb(i)*sin(theta0)-zgb*cos(theta0) )*wi(j) 
0 1 0]; 
elseif delta(k) < 0 
wO = -(awtsqrt (zw~2+4*cdz*E0*zd*delta(k)))/(-2*cdz*E0); 
coeffi = (mw*wO0+cdz*w0"2*E1l+md*delta(k))/wi(j); 
poly1 ((xgb(i)"2+zgb"2) (-2*coeffi*zgb) (coeff1~2-xgb(i)~2)]; 
thetaO = asin(roots(poly1)); 
checki = xgb(i)*cos(theta0(1))+zgb*sin(theta0(1)); 
check2 = xgb(i)*cos(theta0(2))+zgb*sin(theta0(2)); 
if abs(check2-coeffi) < abs(checki-coeffi) 
thetaO = theta0(2); 


else 
thetaO = theta0(1); 
end 
B=[{ (m-zwdot) -(m*xg(i)+zqdot) 0 
-(mwdot+m*xg(i)) iy-mqdot 0) 
0 0 ail 
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A= ((zw-2*cdz*E0*w0) (zq+m) +2*cdz*E1*w0 0 
mw+2*cdz*E1*wO (mq-m*xg(i))-m*zg*w0-2*cdz*E2*w0 


end 


evals = 


end 
end 


(xgb(i)*sin(theta0)-zgb*cos(theta0) )*wi(j) 
1 0); 


eig(A,B); 
degstab(i,j) = max(real(evals)); 


4 Evaluation of the critical speed, ucrit: 


for i=1:length(xg) 
for j=i:(length(u0))-1 


flagi 
flag2 


sign(degstab(i,j)); 
sign(degstab(i,j+1)); 


if flagi “= flag2 


ucrit(i,k) = (u0(j+1)*degstab(i,j)-u0(j)*degstab(i,j+1))/... 
(degstab(i,j)-degstab(i,j+1)); 
end 
end 
end 
end 


4 Revaluation of the nominal angle, ’thetacr’ 
4, at the critical speed, ’ucrit’: 

7, Note that in this case, at a certain ’xg’, 
4 corresponds a specific ’ucrit’. 

7, That’s why we use the same index ’i’: 


w2 = 1556.2363./(nd1.*ucrit.72); 
for i=1:length(xg) 


if delta(1) > 0 
wO = -(zwtsqrt (zw°2-4*cdz*E0*zd*delta(1)))/(2*cdz*E0) ; 
coeff2 = (mw*w0-cdz*w0O~2*E1+md*delta(1i))/w2(i);: 


poly2 = [(xgb(i)"2+zgb°2) (-2*coeff2*zgb) (coeff2°2-xgb(i)~2)]; 
thetaO = asin(roots(poly2)); 
check1 = xgb(i)*cos(theta0(1))+zgb*sin(theta0(1)); 
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check2 = xgb(i)*cos(theta0(2))+zgb*sin(theta0(2)) ; 
if abs(check2-coeff2) < abs(check1-coeff2) 
thetacr(i) = theta0(2); 
else 
thetacr(i) = theta0(1); 
end 
elseif delta(1) < 0 
wO = -(zwtsqrt(zw~2+4*cdz*E0*zd*delta(1)))/(-2*cdz*E0) ; 
coeff2 = (mw*wO+cdz*w0~ 2*E1+md*delta(1))/w2(i):; 
poly2 [(xgb(i)“2+zgb°2) (-2*coeff2*zgb) (coeff2°2-xgb(i)~2)]; 
thetaO = asin(roots(poly2)); 
check1 = xgb(i)*cos(theta0(1))+zgb*sin(theta0(1)); 
check2 = xgb(i)*cos(theta0(2))+zgb*sin(theta0(2)); 
if abs(check2-coeff2) < abs(checki-coeff2) 
thetacr(i) = theta0(2); 
else 
thetacr(i) = theta0(1); 
end 
end 


end 


results = [xg’ ucrit(:,1) wO*ones(length(xg),1) thetacr’] 
save ucritO.dat results -ascii 
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aqQaana 


nm WN 


NOOO fF W DY FE 


PROGRAM HOPF_O.FOR 
EVALUATION OF HOPF BIFURCATION FORMULAS 
ZERO DIVE PLANE ANGLE 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION L,IY,MASS ,MQDOT,MWDOT,ND1,THETA, 
MQ,MW,MD,MDS,MDB,K1,K2, 
ALFA,BETA,GAMA,DELTA, 

EO ,E1,E2,E3,E4, 
DW1,DW2,DW3 ,DW4, 
DQ1,DQ2,D03,D04 

DOUBLE PRECISION M11,M12,M13,M21 ,M22,M23, 
M31,M32,M33, 
N11,N12,N13,N21,N22,N23, 
N31,N32,N33, 
L21,L22,L23,L24,1L31,L32,L33,L34, 
L25,L26,L27,L35,L36,L37, 
L21A,L22A,L23A,L24A,L314A, 
L32A,L33A,L34A 


DIMENSION AC3,3),T(3,3), TINV(3,3) ,FV1(3) ,IV1(3) , YYY(3,3) 
DIMENSION WR(3),WI(3),TSAVE(3,3),TLUD(3,3),IVLUD(3) , SVLUD(3) 
DIMENSION ASAVE(3,3),A1(3,3),A2(3,3) ,XL(25) , BR(25) 

DIMENSION VECO(28),VEC1(25) , VEC2(25) , VEC3(25) , VEC4(25) 

OPEN (20,FILE=’HOPF.RES’ ,STATUS=’ NEW’ ) 


WEIGHT=1556 . 2363 


IY = §61.32 
L = 13.9792 
RHO = 1.94 


WRITE (*,*) ’ ENTER CD’ 
READ (*,*) CD 


CD = 0.5*CD*RHO 
G = 32.2 
XB = 0.0 


WRITE (*,*) ’ ENTER ZG/L’ 
READ (*,*) ZG 

ZG =ZG*L 

ZB =O 

MASS = WEIGHT/G 
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2 


BOY 

wy 8 86="0 
ZQDOT =-6 
ZWDOT =-1 
ZQ = 4 
ZW =-1 
MQDOT =-8 
MWDOT =-5 
MQ =-3 
MW ye 
ZGB=ZG-ZB 
DEFINE THE 
XL( 1)= O 
XL( 2)= 0O 
XL( 3)= 0O 
XL( 4)= 0 
XL( 5)= O 
XL( 6)= 0O 
XL( 7)= O 
XL( 8)= 0 
XL( 9)= 1 
XL(10)= 2 
XL(11)= 3 
XL(12)= 4 
013) 
XL(14)= 10 
XL(15)= 15 
XL(16)= 16 
XL(17)= 17 
XL(18)= 18 
XL(19)= 19 
XL(20)= 20 
XL(21)= 20 
XL(22)= 20 
XL(23)= 20 
XL(24)= 20 
XL(25)= 20 
DO 102 N = 


= WEIGHT 

.5*RHO*L**2 
.3300E-04*0. 
.4629E-02*0 
.5450E-03%*0. 
.3910E-02*0. 
.8000z£-04*0 
.6100E-04*0 
. 1020E-03*0 


0324E-02*0 


LENGTH X AND BREADTH B TERMS FOR THE INTEGRATION 


.0000 
. 1000 
. 2000 
. 3000 
. 4000 
. 8000 
.6000 
. 7000 
0000 
.0000 
.0000 
0000 
. 7143 
.0000 
. 1429 
.0000 
.0000 
.0000 
.0000 
0000 
. 1000 
. 2000 
. 3000 
- 4000 
4167 


1,25 


S5*RHO*L**4 


.5*RHO*L**3 


S*RHO*L**3 
5*RHO*L**2 


.5*RHO*L**5 
.5*RHO*L**4 
.5*RHO*L**4 
-5*RHO*L* *3 
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XL(N) = (L/20.)*XL(N) - L/2. 
102 CONTINUE 


BR( 1)= 0.000 
BR( 2)= 0.485 
BR( 3)= 0.658 
BR( 4)= 0.778 
BR( 5)= 0.871 
BR( 6)= 0.945 
BR( 7)= 1.010 
BR( 8)= 1.060 
BR( 9)= 1.180 
BR(10)= 1.410 
BR(11)= 1.570 
BR(12)= 1.660 
BR(13)= 1.670 
BR(14)= 1.670 
BR(15)= 1.670 
BR(16)= 1.630 
BR(17)= 1.370 
BR(18)= 0.919 
BR(19)= 0.448 
BR(20)= 0.195 
BR(21)= 0.188 
BR(22)= 0.168 
BR(23)= 0.132 
BR(24)= 0.053 
BR(25)= 0.000 


DO 104 K = 1,25 
VECO (K)=BR(K) 
VEC1(K)=XL(K) *BR(K) 
VEC2(K) =XL(K) *XL(K) *BR(K) 
VEC3(K) =XL(K) *XLCK) *XL(K) *BR(K) 
VEC4 (K) =XL(K) *XL(K) *XL(K) *XL(K) *BR(K) 

104 CONTINUE 

CALL TRAP(25,VECO,XL,EO) 

CALL TRAP(25,VECi,XL,E1) 

CALL TRAP(25,VEC2,XL,E2) 

CALL TRAP(25,VEC3,XL,E3) 

CALL TRAP(25,VEC4,XL,E4) 
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Q 


WRITE (*,*) ’ ENTER GAMMA’ 
READ (*,*) EPSILON 
XGMIN=-0.01 

XGMAX=+0.01 

IXG=80 

XGMIN=XGMIN*L 
XGMAX=XGMAX*L 


DO 1 IT=1,IXG 

WRITE (*,3001) IT,IXG 
XG=XGMIN+ (XGMAX-XGMIN) * (IT-1) /(IXG-1) 
XGB=XG-XB 
DV=(MASS-ZWDOT) *( IY-MQDOT) - (MASS*XG+ZQDOT) *(MASS*XG+MWDOT) 
CD6=CD/(3.DO*EPSILON*DV) 
DW1=CD6*( (IY-MQDOT) *(-EO)+(MASS*XG+ZQDOT) *E1) 
DW2=CD6*( (IY-MQDOT) *(3*E1) - (MASS*XG+ZQDOT) *3*E2) 
DW3=CD6* ( (IY-MQDOT) *(-3%*E2) +(MASS*XG+ZQDOT) *3*E3) 
DW4=CD6*( (IY-MQDOT) *(E3) -(MASS*XG+ZQDOT) *E4) 
DQ1=CD6*( (MASS-ZWDOT) *(E1) +(MASS*XG+MWDOT) *(-EO) ) 
DQ2=CD6* ( (MASS-ZWDOT) * (-3*E2) +(MASS*XG+MWDOT) *(3*E1)) 
DQ3=CD6* ( (MASS-ZWDOT ) * (3*E3) +(MASS*XG+MWDOT) * (-3*E2) ) 
DQ4=CD6*( (MASS-ZWDOT) *(-E4)+(MASS*XG+MWDOT) *(E3) ) 
THETAO=ATAN(-XGB/ZGB) 
AAO=(MASS-ZWDOT) * (IY-MQDOT) - (MWDOT+MASS*XG) * (ZQDOT+MASS*XG) 
BBO=(-ZWDOT*MASS-MASS*MW-ZQ*MASS) *XG 

& +(-MASS*MQ+ZWDOT*MQ-ZQDOT*MW 

& -ZQ*MWDOT-MASS*MWDOT-IY*ZW+MQDOT*ZW) 
CCO=-MASS * ZW* XG+MQ* ZW-ZQ*MW-MASS *MW 
CC1=(-MASS *XG+ZWDOT*XG+MASS*XB-ZWDOT*XB)*SIN(THETAO) 


& +(-MASS*ZB-ZWDOT*ZG+ZWDOT*ZB+MASS*ZG) *COS (THETAO) 
DD1=(ZW*XG-ZW*XB) *SIN(THETAO) +(ZW*ZB-ZW*ZG) 
& *COS(THETAO) 


WEI=BB0*CCO/ (AA0*DD1-BBO*CC1) 
UO=DSQRT(1556 .2363/WEI) 
U=U0 


DETERMINE [A] AND [B] COEFFICIENTS 
A11DV=(IY-MQDOT) *ZW+ (MASS*XG+ZQDOT) *MW 


A12DV=(CIY-MQDOT) * (MASS+ZQ) + (MASS *XG+ZQDOT) * (MQ-MASS*XG) 
A13DV=- (MASS*XG+ZQDOT) *WEIGHT 
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A21DV=(MASS-ZWDOT) *MW+ (MASS*XG+MWDOT) *ZW 
A22DV=(MASS-ZWDOT ) * (MQ-MASS*XG)+(MASS*XG+MWDOT) * (MASS+ZQ) 
A23DV=- (MASS-ZWDOT) *WEIGHT 


A11=A11DV/DV 
A12=A12DV/DV 
A13=A13DV/DV 
A21=A21DV/DV 
A22=A22DV/DV 
A23=A23DV/DV 


C11DV=(IY-MQDOT) *MASS*ZG 
C12DV=- (MASS*XG+ZQDOT) *MASS*ZG 
C21DV=- (MASS-ZWDOT) *MASS*ZG 
C22DV=(MASS*XG+MWDOT) *MASS*ZG 


C11=C11DV/DV 
C12=C12DV/DV 
C21=C21DV/DV 
C22=C22DV/DV 


EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 


Ki=-(XGB*SIN(THETAO)-ZGB*COS(THETAO) ) 
K2=-(1./6.)*(ZGB*COS(THETAO)-XGB*SIN(THETAO) ) 


A(1,1)=0.0 
A(1,2)=0.0 
A(1,3)=1.0 
A(2,1)=A13*K1 
A(2,2)=A11*U 
A(2,3)=A12*U 
AC3,1)=A23*K1 
A(3,2)=A21*U 
A(3,3)=A22+U 
DO’ 11 I=1,3 
DO 12 J=1,3 
ASAVE(I,J)=A(I,J) 
i CONTINUE 
11 CONTINUE 
CALL RG(3,3,A,WR,WI,1,YYY,IV1,FV1, IERR) 
CALL DSOMEG(IEV,WR,WI, OMEGA, CHECK) 
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OQ 


14 


13 


16 
17 


OMEGAO=OMEGA 
DO 6° I=1,3 
T(I,1)= YYYC(I,IEV) 
T(I,2)=-YYY(I, IEV+1) 
CONTINUE 
IF (IEV.EQ.1) GO TO 13 
IF (IEV.EQ.2) GO TO 14 
STOP 3004 
DO 6 I=1,3 
Tole sy—11Y(1,1) 
CONTINUE 
GO TO 17 
DO 16 I=1,3 
T(I,3)=YYY(I,3) 
CONTINUE 
CONTINUE 


NORMALIZATION OF THE CRITICAL EIGENVECTOR 


INORM=1 
IF (INORM.NE.O) CALL NORMAL(T) 


INVERT TRANSFORMATION MATRIX 


DO 2 I=1,3 
DO 3 J=1,3 
TINV(I,J)=0.0 
TSAVE(I,J)=T(I,J) 
CONTINUE 
CONTINUE 
CALL DLUD(3,3,TSAVE,3,TLUD,IVLUD) 
DO 4 I=1,3 
IF (IVLUD(I).EQ.0) STOP 3003 
CONTINUE 
CALL DILU(3,3,TLUD,IVLUD,SVLUD) 
DO 8 I=1,3 
DO 9 J=1,3 
TINV(I,J)=TLUD(I,J) 
CONTINUE 
CONTINUE 


CHECK Inv(T)*A*T 
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IMULT=1 

IF (IMULT.EQ.1) CALL MULT(TINV,ASAVE,T,A2) 
IF CIMULT.EQ.0) STOP 

P=A2(3,3) 

PEIG=P 


DEFINITION OF Nij 


Nig =TiNVer 1) 
N21=TINV(2,1) 
N31=TINV(3,1) 
Ni2=TINV G2) 
N22=TINV(2,2) 
N32=TINV(3,2) 
N13=TINV(1,3) 
N23=TINV(2,3) 
N33=TINV(3,3) 


DEFINITION OF Mij 


M11=T(1,1) 
M21=T(2,1) 
M31=T(3,1) 
M12=T(1,2) 
M22=T(2,2) 
M32=T(3,2) 
M13=T(1,3) 
M23=T(2,3) 
M33=T(3,3) 


DEFINITION OF Lij 


L25=C11*M31*M31+C12*M21*M31 
L26=2*C11*M31*M32+C12*(M21*M32+M22+*M31 ) 
L27=C11*M32*M32+C12*M22*M33 
L35=C22*M31*M31+C21*M21¥*M31 
L36=2*C22*M31*M32+C21* (M21*M32+M22*M31 ) 
L37=C22*M32*M32+C21*M33*M22 


DEFINITION OF ALFA, BETA, GAMA 
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D1 =N32*L25 + N33*L35 
D2 =N32*L26 + N33*L36 
D3 =N32*L27 + N33*L37 


D1ii=-P 
D1i2=OMEGAO 
D21=-2*0MEGAO 
D22=-P 
D23=2*O0MEGAO 
D32=-OMEGAO 
D33=-P 


BETA=(D2-D21*D1/D11-D23*D3/D33) /(D22-D21*D12/D11-D23*D32/D33) 
ALFA=(D1-D12*BETA)/D11 

GAMA=(D3-D32*BETA)/D33 
L21A=2*C11*ALFA*M31*M33+C12*ALFA*(M21%*M33+M23*M31 ) 
L22A=2*C11*ALFA*M32*M33 + 2*C11*BETA*M31*M33 

+ C12*ALFA*(M22*M33+M32*M23) 

+ (C12*BETA*(M21*M33+M23*M31) 
L23A=2*C11*GAMA*M31*M33+2%*C11*BETA*M32*M33 

+ C12*GAMA*(M21*M33+M23%M31) 

+ (C12*BETA*(M22*M33+M23*M32) 
L24A=2*C11*GAMA*M32*M33+C12*GAMA*(M22*M33+M23*M32) 
L31A=2*C22*ALFA*M31*M33+C21*ALFA*(M21*M334+M23%M31 ) 
L32A=2*C22*ALFA*M32*M33+2*C22*BETA*M31*M33 

+ (C21*ALFA* (M22*M33+M32*M23) 

+ C21*BETA*(M21*M33+M23*M31 ) 
L33A=2*C22*GAMA*M31*M33+2*C22*BETA*M32*M33 

+ (C21*GAMA*(M21*M33+M23*M31) 

+ C21*BETA* (M22*M33+M23*M32) 
L34A=2*C22*GAMA*M32*M334+C21*GAMA* (M22*M33+M23*M32) 
L21=L21A+A13*K2*M11%**3+DW1*M21%*3+DW2*M31*M21%«2 

+ DW3*M21*M31**2+DW4#M31%%3 
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L22=L22A+3*A13*K2*M12*M11**2+3*DW1*M22*M21 **2 


+ DW2* (2*M21*M22*M31+M32*M21 **2 ) 

+ DW3* (2*M21*M31*M32+M22*M31**2) 

+ 3*DW4*M32*M31**2 
L23=L23A+3*A13*K2*M11*M12**2+3%DW1*M21*M22*%2 

+ DW2* (M31*M22* *2+2*M21*M22*M32) 

+ DW3* (M21*M32**2+2%*M22*M31*M32) 

+ 3*DW4*M31*M32**2 
L24=L24A+A13*K2*M12**3+DW1*M22**3+DW2*M32*M22*%2 

+ DW3*M22*M32**2+DW4*M32**3 
L31=L31A+A23*K2*M11**3+DQ1*M21**3+DQ2*M31*M21*%2 

+ DQ3*M21*M31%**2+DQ4*M31**3 
L32=L32A+3*A23*K2*M12*M11%**2+3*DQ1*#M22*M21%**2 

+ DQ2* (2*M21*M22*M31+M32*M21**2 ) 

“* DQ3* (2*M21*M31+*M32+M22*M31%**2) 

+ 3*DQ4*M32*M31**2 
L33=L33A+3*A23*K2*M11*M12**2+3*DQ1*M21*M22**2 

+ DQ2* (M31*M22**2+2*M21*M22*M32) 

+ DQ3* (M21*M32**2+2*M22*M31*M32) 

+ 3*DQ4*M31*M32**2 
L34=L34A+A23 *K2*M12**3+DQ1*M22**3+DQ2*M32*M22**2 

+ DQ3*M22*M32**2+DQ4*M32*%3 


R1i1=N12*L21+N13*L31 
R12=N12*L22+N13*L32 
R13=N12*L23+N13*L33 
R1i4=N12*L24+N13*L34 
R21=N22*L21+N23*L31 
R22=N22*L22+N23*L32 
R23=N22*L23+N23*L33 
R24=N22*L24+N23*L34 


EVALUATE DALPHA AND DOMEGA 


UINC=0.001 
UR =U+UINC 
UL =U-UINC 
U  =UR 


A(1,1) 
A(1,2) 
cG.s) 


0.0 
0.0 
1.0 
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AC 2,1)=A13*K1 
A(2,2)=A11*U 
A(2,3)=A12*U 
A(3,1)=A23*K1 
A(3,2)=A21*U 
AC3,3)=A22*U 


CALL RG(3,3,A,WR,WI,0,YYY,IV1,FV1, TERR) 
CALL DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEO0S 

OMEGR=FREQ 


=UL 


A(1,1)=0.0 
A(1,2)=0.0 
AC1,3)=1.0 
A(2,1)=A13*K1 
AC2,2)=A11*U 
A(2,3)=A12*U 
A(3,1)=A23*K1i 
AC3,2)=A21*U 
A(3,3)=A22*U 


CAEL RG(3,3,A,WR,WL,O,YYY,IV1,FV1,IERR) 
CALL DSTABL(DEOS ,WR,WI,FREQ) 

ALPHL=DE05S 

OMEGL=FREQ 


DALPHA=(ALPHR-ALPHL) / (UR-UL) 
DOMEGA=(OMEGR-OMEGL) / (UR-UL) 


EVALUATION OF HOPF BIFURCATION COEFFICIENTS 


COEF 1=3 .0*R11+R13+R22+3 .0*R24 

COEF2=3 .0*R21+R23-R12-3.0*R14 

AMU2 =-COEF1/(8.0*DALPHA) 

BETA2=0 .25*COEF1 

TAU2 =-(COEF2-DOMEGA*COEF1/DALPHA)/(8.0*0MEGAO) 
PER =2.0*3.1415927/0MEGAO 

PER =PER*U/L 

WRITE (20,2001) XG/L,U,COEF1,DALPHA,OMEGAO,PEIG 
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1 CONTINUE 
STOP 

1001 FORMAT (’ ENTER NUMBER OF DATA LINES’ ) 
1002 FORMAT (’ ENTER UO, ZG, AND DSAT’) 
1003 FORMAT (’ ENTER BOW PLANE TO STERN PLANE RATIO’) 
1004 FORMAT (’ ENTER ZG’) 
2001 FORMAT (6E14.5) 
4001 FORMAT (3F15.5) 
3001 FORMAT (215) 
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SUBROUTINE DSOMEG(IJK,WR,WI,OMEGA ,CHECK) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(3),WI(3) 
CHECK=-1.0D+25 
DO 1 I=1,3 
IF (WR(I).LT.CHECK) GO TO 1 
CHECK=WR(I) 
1S 
1 CONTINUE 
OMEGA=DABS(WICIJ)) 
IF (WI(IJ).GT.0.DO0) IJK=IJ 
IF (WI(IJ).LT.0.DO) IJK=IJ-1 
RETURN 


SUBROUTINE DSTABL(DEOS,WR, WI, OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(3),WI(3) 
DEOS=-1.0D+20 
DO 1 I=1,3 
IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR(I) 
IJzI 
1 CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS (OMEGA) 
RETURN 


SUBROUTINE NORMAL(T) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
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DIMENSION T(3,3),TNOR(3,3) 
CFAC=T(1,1)#*2+T(1,2)**2 
IF (DABS(CFAC).LE.(1.D-10)) STOP 4001 
TNOR(1,1)=1.D0 
MmiOR(2,1)=(T(1 ,1)*1T(2,1)+T(2,2)*T(1,2))/CFAC 
TNOR(3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2))/CFAC 
TNOR(1,2)=0.D0 
TNOR(2,2)=(T(2,2)*T(1,1)-T(2,1)*T(1,2))/CFAC 
TNOR(3,2)=(T(3,2)*T(1,1)-T(3,1)*T(1,2))/CFAC 
DO 1 I=1,3 

Dole2 J=1.2 

TOT , J) =LNORC IO 
CONTINUE 
1 CONTINUE 

RETURN 


SUBROUTINE MULT(TINV,A,T,A2) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION TINV(3,3),A(3,3),T(3,3) ,A1(3,3) ,A2(3,3) 
Dod I=1,3 
DOM2 J=1,3 
A1(I,J)=0.DO 
A2(1I,J)=0.DO 
2 CONTINUE 
1 CONTINUE 
DO 3 I=1,3 
DO 4 J=1,3 
DO & K=1,3 
A1(I,J)=ACI,K)*T(K,J)+A1(1I,J) 
5 CONTINUE 
CONTINUE 
3 CONTINUE 
DO 6 I=1,3 
DO 7 J=1,3 
DO 8 K=1,3 
A2(I,J)=TINV(I,K)*A1(K,J)+A2(1I,J) 
8 CONTINUE 
CONTINUE 
6 CONTINUE 
DO 11 I=1,3 
GC WRITE (*,101) (ACI,J),J=1,3) 


iD 
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CONTINUE 
DO 12 I=1,3 
WRITE (*,101) (T(I,J),J=1,3) 
CONTINUE 
DO 10 I=1,3 
WRITE (*,101) (A2(I,J),J=1,3) 
CONTINUE 
WRITE (*,101) A2(1,1) 
RETURN 
FORMAT (3E15.5) 
END 


SUBROUTINE TRAP(N,A,B,OUT) 
NUMERICAL INTEGRATION ROUTINE USING THE TRAPEZOIDAL RULE 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION A(1),B(1) 

Ni=N-1 

OUT=0.0 

DO 1 I=1,N1 
OUT1=0.5*(ACI)+ACI+1))*(BCI+1)-B(I)) 
OUT =OUT+OUT1 

CONTINUE 

RETURN 

END 
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PROGRAM HOPF_DELTA.FOR 
EVALUATION OF HOPF BIFURCATION FORMULAS 
STERN PLANE ANGLE IS NON ZERO 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION L,IY,MASS,MQDOT ,MWDOT,ND1, THETA, 
MQ,MW,MD,ZD,MDS ,MDB,K1,K2, 
ALFA,BETA,GAMA,DELTA, 

EO ,E1,E2,E3,E4, 
DW1,DW2,DW3 ,DW4, 
DQ1,DQ2,DQ3 ,DQ4 

DOUBLE PRECISION M11 ,Mi2,M13,M21,M22,M23, 

M31 ,M32,M33, 
N11,N12,N13 ,N21 N22 ,N23; 
NSi,N3S2,N33, 
L21,L22,L23,L24,L31,L32,L33 ,L34, 
L25,L26 ,L27 ,L35,L36,L37, 

L214 ,L22A,L23A,L244,L314, 
L32A4,L334,L344 


DIMENSION A(3,3),T(3,3),TINV(3,3) ,FV1(3) ,IV1(3) , YYY(3,3) 
DIMENSION WR(3),WI(3) ,TSAVE(3,3) ,TLUD(3,3),IVLUD(3) , SVLUD(3) 
DIMENSION ASAVE(3,3),A1(3,3),A2(3,3) ,XL(25) ,BR(25) 

DIMENSION VECO(25),VEC1(25),VEC2(25) , VEC3(25) , VEC4(25) 


OPEN (10,FILE=’HOPF .DAT’ ,STATUS=’OLD’ ) 
OPEN (20,FILE=’HOPF.RES’ ,STATUS=’ NEW’ ) 


WEIGHT=1556 . 2363 


ey = 661.32 
L = 13.9792 
RHO = 1.94 


WRITE (*,*) ’ ENTER CD’ 
READ (*,*) CD 


CD = 0.5*CD*RHO 
G Seo. 2 
XB = 0.0 


WRITE (*,*) ’ ENTER ZG/L’ 
READ (*,*) ZG 

ZG =ZG*L 

ZB = 0.0 
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MASS = WEIGHT/G 
BOY = WEIGHT 
ND1 = 0.5*RHO*L**2 


ZQDOT =-6.3300E-04%*0.5*RHO*L**4 
ZWDOT =-1.4529E-02*0.5*RHO*L**3 
ZQ = 7.5460E-03*0 .6*RHO*L**3 
ZW =-1.3910E-02*0.5*RHO*L**2 
MQDOT =-8.8000E-04*0.5*RHO*L**5& 
MWDOT =-&.6100E-04*0 .5*RHO*L**4 


MQ =-3.7020E-03*0.5*RHO*L**4 
MW = 1.0324E-02*0.5*RHO*L**3 
ZGB=ZG-ZB 


DEFINE THE LENGTH X AND BREADTH B TERMS FOR THE INTEGRATION 


XL( 1)= 0.0000 
XL( 2)= 0.1000 
XL( 3)= 0.2000 
XL( 4)= 0.3000 
XL( 5)= 0.4000 
XL( 6)= 0.5000 
XL( 7)= 0.6000 
XL( 8)= 0.7000 
XL( 9)= 1.0000 
XL(10)= 2.0000 
XL(i1)= 3.0000 
XL(12)= 4.0000 
XL(13)= 7.7143 


XL(14)= 10.0000 
XL(16)= 15.1429 
XL(16)= 16.0000 
XL(17)= 17.0000 
XL(18)= 18.0000 
XL(19)= 19.0000 
XL(20)= 20.0000 
XL(21)= 20.1000 
XL(22)= 20.2000 
XL(23)= 20.3000 
XL(24)= 20.4000 
XL(25)= 20.4167 
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DO 102 N sei,25 
XL(N) = (L/20.)*XL(N) - L/2. 
102 CONTINUE 


BR( 1)= 0.000 
BR( 2)= 0.485 
BR( 3)= 0.658 
BRC 4)= 0.778 
BRC 5)= 0.871 
BR( 6)= 0.945 
BR( 7)= 1.010 
BR( 8)= 1.060 
BRC 9)= 1.180 
BR(10)= 1.410 
BR(11)= 1.570 
BR(12)= 1.660 
BR(13)= 1.670 
BR(14)= 1.670 
BR(15)= 1.670 
BR(16)= 1.630 
BR(17)= 1.370 
BR(18)= 0.919 
BR(19)= 0.448 
BR(20)= 0.195 
BR(21)= 0.188 
BR(22)= 0.168 
BR(23)= 0.132 
BR(24)= 0.053 
BR(25)= 0.000 


DO 104 K = 1,25 
VECO (K)=BR(K) 
VEC1(K)=XL(K) *BR(K) 
VEC2(K)=XL(K) *XL(K) *BR(K) 
VEC3 (K) =XL(K) *XL(K) *XL(K) *BR(K) 
VEC4(K)=XL(K) *XL(K) *XL(K) *XL(K) *BR(K) 
104 CONTINUE 
CALL TRAP(25,VECO,XL,EO) 
CALL TRAP(25,VEC1,XL,E1) 
CALL TRAP(25,VEC2,XL,E2) 
CALL TRAP(25,VEC3,XL,E3) 
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CALL TRAP(25,VEC4,XL,E4) 
WRITE (*,*) ’ ENTER GAMMA’ 
READ (*,*) EPSILON 


DO 1 IT=1;1%E 
READ (10,*) XG,U,WO,THETAO 
WRITE (*,3001) IT,IXG 
XG=XG*L 
XGB=XG-XB 
DV=(MASS-ZWDOT) * (IY-MQDOT) - (MASS*XG+ZQDOT) * (MASS*XG+MWDOT) 
CONST2=CD/ (3.0*DV*EPSILON) 
TDW1=CONST2*( (IY-MQDOT)»*(-EO)+(MASS*XG+ZQDOT) *(E1) ) 
TDW2=CONST2* ( (IY-MQDOT) *(3.0*E1)+(MASS*XG+ZQDOT) *(-3.0*E2) ) 
TDW3=CONST2*( (LY-MQDOT) *(-3.0*E2)+(MASS*XG+ZQDOT)*(3.0*E3)) 
TDW4=CONST2*( (IY-MQDOT) *(E3)+(MASS*XG+ZQDOT) *(-E4) ) 
TDQ1=CONST2*( (MASS-ZWDOT )* (E1) + (MASS *XG+MWDOT) *(-E0) ) 
TDQ2=CONST2* ( (MASS-ZWDOT) *(-3.0*E2)+(MASS*XG+MWDOT) *(3.0*E1) ) 
TDQ3=CONST2* ( (MASS-ZWDOT) *(3.0*E3) +(MASS*XG+MWDOT) *(-3.0*E2)) 
TDQ4=CONST2*( (MASS-ZWDOT) *(-E4)+(MASS*XG+MWDOT) *(E3) ) 


CONST=TANH(WO/EPSILON) *CD/DV 

DW1=CONST*( CIY-MQDOT) *(-E0) +(MASS*XG+ZQDOT) *E1) 
DW2=CONST*( (IY-MQDOT)*(2*E1) -(MASS*XG+ZQDOT) *2*E2) 
DW3=CONST*( (CIY-MQDOT)*(-E2) +(MASS*XG+ZQDOT) *E3) 
DQi=CONST*( (MASS-ZWDOT)*(E1)  +(MASS*XG+MWDOT)*(-EO) ) 
DQ2=CONST*( (MASS-ZWDOT) *(-2*E2)+(MASS*XG+MWDOT) *2*E1 ) 
DQ3=CONST*( (MASS-ZWDOT)*(E3)  +(MASS*XG+MWDOT) *(-E2) ) 


DETERMINE [A] AND [B] COEFFICIENTS 


A11DV=(IY-MQDOT) * (ZW-2*CD*E0*WO) 
& +(MASS*XG+ZQDOT) * (MW+2*CD*E1*WO) 
A12DV=(IY-MQDOT ) *(MASS+ZQ+2*CD*E1*WO) +(MASS*XG+ZQDOT) 
& *(MQ-MASS*XG-MASS*ZG*WO-2*CD*E2*WO) 
A13DV=-(MASS*XG+ZQDOT) *WEIGHT 
A21DV=(MASS-ZWDOT) * (MW+2*CD*E1*WO) 


& +(MASS*XG+MWDOT) *(ZW-2*CD*EO*WO ) 
A22DV=(MASS-ZWDOT)* (MQ-MASS*XG-MASS*ZG*WO-2*CD*E2*WO) 
& +(MASS*XG+MWDOT) * (MASS+ZQ+2*CD*E1*WO) 


A23DV=-(MASS-ZWDOT) *WEIGHT 
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A11=A11DV/DV 
A12=A12DV/DV 
A13=A13DV/DV 
A21=A21DV/DV 
A22=A22DV/DV 
A23=A23DV/DV 


C11DV=(IY-MQDOT) *MASS*ZG 
C12DV=-(MASS*XG+ZQDOT) *MASS*ZG 
C21DV=-(MASS-ZWDOT) *MASS*ZG 
C22DV=(MASS*XG+MWDOT ) *MASS*ZG 


C11=C11DV/DV 
C12=C12DV/DV 
C21=C21DV/DV 
C22=C22DV/DV 


EVALUATE TRANSFORMATION MATRIX OF EIGENVECTORS 


Ki=-(XGB*SIN(THETAO) -ZGB*COS (THETAO) ) 
K2=-(1./6.)*(ZGB*COS(THETAO) -XGB*SIN(THETAO) ) 


A(1,1)=0.0 
A(1,2)=0.0 
A(1,3)=1.0 
AC2,1)=A13*K1 
A(2,2)=A11*U 
AC2,3)=A12#U 
A(3,1)=A23*K1 
A(3,2)=A21+U 
A(3,3)=A22+*U 
Berit I=1,3 
DO 12 J=1,3 
ASAVE(I,J)=ACI,J) 
12 CONTINUE 
11 CONTINUE 
CALL RG(3,3,A,WR,W1I,1,YYY,IV1,FVi, IERR) 
CALL DSOMEG(IEV,WR,WI,OMEGA, CHECK) 
OMEGAO=O0MEGA 
DO § I=1,3 
T(I,1)= YYY(1,IEV) 
T(I,2)=-YYY(I, IEV+1) 
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CONTINUE 
IF (IEV.EQ.1) GO TO 13 
IF (IEV.EQ.2) GO TO 14 
STOP 3004 
DO 6 I=1,3 
TCI,3)=YYY(I,1) 
CONTINUE 
GO TO 17 
DO 16 I=1,3 
ECL sari (Los) 
CONTINUE 
CONTINUE 


NORMALIZATION OF THE CRITICAL EIGENVECTOR 


INORM=1 
IF (INORM.NE.O) CALL NORMAL(T) 


INVERT TRANSFORMATION MATRIX 


DO 2 I=1,3 
DO 3 J=1,3 
TINV(I,J)=0.0 
TSAVE(I,J)=T(I,J) 
CONTINUE 
CONTINUE 
CALL DLUD(3,3,TSAVE,3,TLUD,IVLUD) 
DO 4 I=1,3 
IF (IVLUD(I).EQ.0) STOP 3003 
CONTINUE 
CALL DILU(3,3,TLUD,IVLUD,SVLUD) 
DO 8 I=1,3 
DO 9 J=1,3 
TINV(I,J)=TLUD(I, J) 
CONTINUE 
CONTINUE 


CHECK Inv(T)*A*T 
IMULT=1 


IF (IMULT.EQ.1) CALL MULTC(TINV,ASAVE,T,A2) 
IF CIMULT.EQ.0) STOP 
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P=A2(3,3) 
PEIG=P 


DEFINITION OF Nij 


N11=TINV(1,1) 
N21=TINV(2,1) 
N31=TINV(3,1) 
N12=TINV(1,2) 
N22=TINV(2,2) 
N32=TINV(3,2) 
N13=TINV(1,3) 
N23=TINV(2,3) 
N33=TINV(3,3) 


DEFINITION OF Mij 


M11=T(1,1) 
M21=T(2,1) 
M31=T(3,1) 
M12=T(1,2) 
M22=T (2,2) 
M32=T(3,2) 
M13=T(1,3) 
M23=T(2,3) 
M33=T(3,3) 


DEFINITION OF Lij 


L25=C11*M31*M31+C12*M21*M31 
L26=2*C11*M31*M32+C12*(M21*M32+M22*M31 ) 
L27=C11*M32*M32+C12*M22*M33 
L35=C22*M31*M31+C21*M21*M31 
L36=2*C22*M31*M32+C21*(M21*M32+M22*M31) 
L37=C22*M32*M32+C21*M33*M22 


DEFINITION OF ALFA, BETA, GAMA 
D1 =N32*L25 + N33*L35 


D2 =N32*L26 + N33*L36 
D3 =N32*L27 + N33*L37 
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Dift==P 
D12=OMEGAO 
D21=-2*0MEGAO 
D22=-P 
D23=2*0MEGAO 
D32=-OMEGAO 
D33=-P 


BETA=(D2-D21*D1/" ** -D23*D3/D33) / (D22-D21*D12/D11-D23*D32/D33) 
ALFA=(D1-D12*BETA)/D11 
GAMA=(D3-D32*BETA)/D33 


L21A=2*(C11+DW3) *ALFA*M31*M33+(C12+DW2) *ALFA*(M21*M33+M23*M31) 
+2*DW1*M21*M23*ALFA 


L22A=2*(C11+DW3) *ALFA*M32*M33 + 2*(C11+DW3) *BETA*M31%*M33 
+ (C12+DW2) *ALFA* (M22*M32+M32*M23) 
+ (C12+DW2)*BETA*(M21*M33+M23¥*M31) 
+ 2*DW1*(M21*M23*BETA+M22*M23*ALFA) 


L23A=2* (C11+DW3) *GAMA*M31*M33+2*(C11+DW3) *BETA*M32*M33 
+ (C12+DW2)*GAMA*(M21*M33+M23*M31) 
+ (C12+DW2)*BETA* (M22*M33+M23*M32) 
+ 2*DW1*(M21*M23*GAMA+M22*M23*BETA) 


L24A=2* (C11+DW3 ) *GAMA*M32*M33+(C12+DW2) *GAMA* (M22*M33+M23*M32) 
+2*DW1*M22*M23*GAMA 


L31A=2* (C22+DQ3) * ALFA*M31*M33+(C21+DQ2) *ALFA* (M21*M33+M23¥*M31 ) 
+2*DQ1*M21*M23*ALFA 


L32A=2% (C22+DQ3) *ALFA*M32*M33 + 2*(C22+DQ3)*BETA*M31*M33 
+ (C21+DQ2) *ALFA* (M22*M33+M32*M23) 
+ (C21+DQ2)*BETA*(M21*M33+M23%M31) 
+ 2*DQ1*(M21*M23*BETA+M22*M23*ALFA) 


L33A=2* (C22+DQ3) *GAMA*M31*#M33+2* (C22+DQ3) *BETA*M32*M33 
+ (C21+DQ2) *GAMA* (M21*M33+M23*M31) 
+ (C21+DQ2) *BETA* (M22*M33+M23*M32) 
+ 2*DQ1*(M21*M23*GAMA+M22*M23*BETA) 


L34A=2* (C22+DQ3) *GAMA*M32*M33+ (C21+DQ2) *GAMA* (M22*M33+M23*M32) 


Q 


& 


+2*DQ1*M22*M23*GAMA 


L21=L21A+A13*K2*M11%*3 
L22=L22A+3*A13*K2*M12*M11**2 
L23=L23A+3*A13*K2*M11*M12**2 
L24=L24A+A13*K24*M12%*3 
L31=L31A+A23*K2¥*M11**3 
L32=L32A+3*A23*K2*M12*M11**2 
L33=L33A+3*A23*K2*M11*M12**2 
L34=L34A+A23*K2*M12%%3 


L21=L21+TDW1*M21**3+TDW2*M31*M21**2+TDW3*M21*M31**2+TDW4*M31**3 
L22=L22+TDW1*3 .0*M22*M21%**2+TDW2* (2 .0#*M21*M31*M22+M324*M21%**2) 
+TDW3*(2.0#M314*M32*M214+M22*M31%**2)+TDW4*3 .0*M32*M31**2 
L23=L23+TDW1*3 .0*M21*M22**2+TDW2*(M31*M22**24+2 .0*M21*M22*M32 ) 
+TDW3* (M21*M32**24+2 .0*M31*M32*M22)+TDW4*3 .0*M31*M32**2 
L24=L244+TDW1*M22**34+TDW2*M32*M22**2+TDW3¥*M22*M32**2+TDW4*M32**3 
L31=L31+TDQ1*M21**3+TDQ2*M31*M21**2+TDQ3*M21*M31**2+TDQ4*M31**3 
L32=L32+TDQ1*3 .0*M22*M21**2+TDQ2*(2.0*M21*M31*M22+M32*M21%*2) 
+TDQ3%*(2.0%*M31*M32*M21+M22*M31¥**2)+TDQ4*3 .0*M324*M31%**2 
L33=L33+TDQ1*3 .0*M21*M22**2+TDQ2*(M31*M22**24+2 .0*M21*M22*M32 ) 
+TDQ3*(M21*M32**24+2 .0*M31*M32*M22 )+TDQ4*3 .0*M31*M32**2 
L34=L34+TDQ1*M22**3+TDQ2*M32*M22**2+TDQ3*M22*M32**2+TDQ4*M32**3 


R11=N12*L21+N13*L31 
R12=N12*L22+N13*L32 
R13=N12*L23+N13*L33 
R14=N12*L244+N13*L34 
R21=N22*L21+N23*L31 
R22=N22*L22+N23*L32 
R23=N22*L23+N23*L33 
R24=N22*L24+N23*L34 


EVALUATE DALPHA AND DOMEGA 
UVINC=0 .001 

UR =U+UINC 

UL =U-UINC 

U =UR 

A(1,1)=0.0 

A(1,2)=0.0 
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A(1,3)=1.0 

AC(2,1)=A13*K1 
AC 2,2)=A11*U 
AC2,3)=A12*U 
AC3,1)=A23*K1 
A(3,2)=A214U 
AC3,3)=A224*U 


CALL RG(3,3,A,WR,WI,0,YYY,IV1,FV1,IERR) 
CALL DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEOS 

OMEGR=FREQ 


U=UL 


A(1,1)=0.0 
A(1,2)=0.0 
A(1,3)=1.0 
A(2,1)=A13*K1i 
AC2,2)=A11*U 
AC 2,3) =A124U 
AC3,1)=A23*K1 
A(3,2)=A21*U 
AC(3,3)=A22*U 


CALL RG(3,3,A,WR,WI,0,YYY,IV1,FV1, IERR) 
CALL DSTABL(DEOS ,WR,WI, FREQ) 

ALPHL=DEOS 

OMEGL=FREQ 


DALPHA=(ALPHR-ALPHL) / (UR-UL) 
DOMEGA=(OMEGR-OMEGL) / (UR-UL) 


EVALUATION OF HOPF BIFURCATION COEFFICIENTS 


COEF 1=3 .0*R11+R13+R22+3 .0*R24 

COEF 2=3 .0*R21+R23-R12-3.0*R14 

AMU2 =-COEF1/(8.0*DALPHA) 

BETA2=0.25*COEF1 

TAU2 =-(COEF2-DOMEGA*COEF1/DALPHA)/(8.0*0MEGAO) 
PER =2.0*3.1415927/0MEGAO 

PER =PER*U/L 
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WRITE (20,2001) XG/L,U,COEF1,DALPHA,OMEGAO,PEIG 
1 CONTINUE 
STOP 
1001 FORMAT (’ ENTER NUMBER OF DATA LINES’) 
1002 FORMAT (’ ENTER UO, ZG, AND DSAT’) 
1003 FORMAT (’ ENTER BOW PLANE TO STERN PLANE RATIO’) 
1004 FORMAT (’ ENTER ZG’) 
2001 FORMAT (6E14.5) 
4001 FORMAT (3F15.5) 
3001 FORMAT (215) 


END 
CHSStSSattSSSSSSrSSesSSSSSSSifsSis2ssssqqsasssssessessrssrrss===z==z= 
SUBROUTINE DSOMEG(IJK,WR,WI,OMEGA,CHECK) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(3) ,WI(3) 
CHECK=-1.0D+25 
DO 1 J=1,3 
IF (WR(I).LT.CHECK) GO TO i 
CHECK=WR(I) 
IJ=I 
1 CONTINUE 
OMEGA=DABS(WI(IJ)) 
IF (WI(IJ).GT.0.DO) IJK=IJ 
IF (WI(IJ).LT.0.DO0) IJK=IJ-1 
RETURN 
END 
Cees SSSSSSSSSesseS2SSSSSSeSSeSseesssssssqsssseesssqssssssess2222=22= 
SUBROUTINE DSTABL(DEOS,WR,WI,OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(3) ,WI(3) 
DEOS=-1 .0D+20 
DO 1 I=1,3 
IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR(I) 
IJ=I 
1 CONTINUE 
OMEGA=WICIJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 
CeSSSSSSS2SSSSSSSSSSSSSSSSSSSSSSSSS2SSSSSSSSSS82SSSSSSSSS 5 2ES=S=S=S== 


SUBROUTINE NORMAL(T) 
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IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION T(3,3),TNOR(3,3) 
CFAC=T(1,1)**2+T(1,2) **2 
IF (DABS(CFAC).LE.(1.D-10)) STOP 4001 
TNOR(1,1)=1.D0 
TNOR(2,1)=(T(1,1)*T(2,1)+T(2,2)*T(1,2))/CFAC 
TNOR(3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2))/CFAC 
TNOR(1,2)=0.D0 
TNOR(2,2)=(T(2,2)*T(1,1)-T(2,1)*T(1,2))/CFAC 
TNOR(3 ,2)=(T(3,2)*T(1,1)-T(3,1)*T(1,2))/CFAC 
DOM I=t2 
DOr 299=1 ,2 
T(I,J)=TNORCI, J) 

2 CONTINUE 

1 CONTINUE 
RETURN 


SUBROUTINE MULT(TINV,A,T,A2) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION TINV(3,3),A(3,3),T(3,3) ,A1(3,3) ,A2(3, 3) 
DO 1 I=1,3 
DO 2 J=1,3 
A1(I,J)=0.D0 
A2(I,J)=0.DO 
CONTINUE 
1 CONTINUE 
DO 3 I-13 
DO 4 J=1,3 
DO § K=1,3 
A1(I,J)=AC(I,K)*T(K,J)+A1(I,J) 
5 CONTINUE 
CONTINUE 
3 CONTINUE 
DO 6 I=1,3 
Dee? J=1;3 
DO 8 K=1,3 
A2(I,J)=TINV(I,K)*A1(K,J)+A2(1I,J) 
8 CONTINUE 
CONTINUE 
6 CONTINUE 
pO 11 I=1,3 


He 


=~] 
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WRITE (*,101) (AC(I,J),J=1,3) 
CONTINUE 
Pur 12 I=1,3 

WRITE (*,101) (T(I,J) ,J=1,3) 
CONTINUE 
Perio I=1,3 

WRITE (*,101) (A2(I,J),J=1,3) 

CONTINUE 
WRITE (*,101) A2(1,1) 
RETURN 
FORMAT (3E15.5) 
END 


SUBROUTINE TRAP(N,A,B,OUT) 
NUMERICAL INTEGRATION ROUTINE USING THE TRAPEZOIDAL RULE 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION A(i),B(1) 

Ni=N-1 

OUT=0.0 

Doma I=4, N41 
OUT1=0.5*( ACI) +ACI+1))*(BCI+1)-B(I)) 
OUT =OUT+0UT1 

CONTINUE 

RETURN 

END 
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PROGRAM SIM.FOR 


SIMULATION PROGRAM USING A FOURTH ORDER ADAMS-BASHFORTH 
INTEGRATION SCHEME 


REAL L,MW,MASS,IY,MQDOT,MWDOT,MQ, 
U,TIME,ND1,ND2,ND3,ND4, 
F1i(4) ,F2(4) ,F3(4) 


DIMENSION XL(25),BR(25S) , VEC1(25) , VEC2(25) 


OPEN (10,FILE=’SIM.DAT’ ,STATUS=’ OLD’ ) 
OPEN (20,FILE=’SIM.RES’ ,STATUS=’ NEW’ ) 


ENTER SPEEDS AND TIME DATA 


READ(10,*) U,TSIM,DELT,IPRNT, 
ZG ,X4G, THETA ,W,Q 


GEOMETRIC PROPERTIES AND HYDRODYNAMIC COEFFICIENTS 


PI =4,.0*ATAN(1.0) 
THETA =THETA*PI/180 


RHO =1.94 
CDZ =0.50 
L =13.9792 


ND1 =0.5*RHO*L**2 

ND2 =0.5*RHO*L**3 

ND3  =0.5*RHO*L**4 

ND4 =0.5*RHO*L**5 
WEIGHT=1556 .2363/ (ND1*U**2) 
MASS =1586.2363/(32.2*ND2) 
IY =661.32/ND4 

ZQDOT=-6 .3300E-04 
ZWDOT=-1.4529E-02 

ZQ = 7.5450E-03 

ZW = =-1.3910E-02 

MQDOT=-8 .8000E-04 
MWDOT=-5.6100E-04 

MQ =-3.7020E-03 
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MW = 1.0324E-02 


ZB =0.0 
XB =0.0 
ZGB =ZG-ZB 
XGB =XG-XB 


DETERMINE [A] AND [B] COEFFICIENTS 


DV =(MASS~ZWDOT) * CLY-MQDOT) ~ (MASS*XG+ZQDOT) * (MASS*XG+MWDOT) 
A11DV=CIY~-MQDOT) *ZW+(MASS*XG+ZQDOT) *MW 
A12DV=(ILY-MQDOT) * (MASS+ZQ) +(MASS*XG+ZQDOT) * (MQ-MASS*XG) 
A13DV=~ (MASS*XG+ZQDOT) *WEIGHT 
A21DV=(MASS-ZWDOT) *MW+ (MASS*XG+MWDOT) *ZW 
A22DV=(MASS-ZWDOT) * (MQ-MASS*XG) + (MASS*XG+MWDOT) *(MASS+ZQ) 
A23DV=~ (MASS~ZWDOT ) *WEIGHT 


A11=A11DV/DV 
A12=A12DV/DV 
A13=A13DV/DV 
A21=A21DV/DV 
A22=A22DV/DV 
A23=A23DV/DV 


DEFINE THE 


XL( 1)= 
Kine 2) 
XL( 3)= 
XL( 4)= 
XL( 5)= 
XL( 6)= 


XL( 8)= 
XL( 9)= 
XL(10)= 
XL(11)= 
Met) = 
XL(13)= 


XL(16)= 16 


LENGTH X AND BREADTH B TERMS FOR THE INTEGRATION 


0.0000 
0.1000 
0.2000 
0.3000 
0.4000 
0.5000 

XL( 7)= 0. 
0 
1 
2 
3 
4 
7 


6000 


. 7000 
.0000 
.0000 
.0000 
.0000 
.7143 
XL(14)= 10. 
XL(15)= 15. 
.0000 


0000 
1429 
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XL(17)= 17.0000 
XL(18)= 18.0000 
XL(19)= 19.0000 
XL(20)= 20.0000 
XL(21)= 20.1000 
XL(22)= 20.2000 
XL(23)= 20.3000 
XL(24)= 20.4000 
XL(25)= 20.4167 


DO 102 N = 1,25 
XL(N) = (L/20)*XL(N) - L/2 
102 CONTINUE 


BRC 1)= 0.000 
BR( 2)= 0.485 
BR( 3)= 0.658 
BR( 4)= 0.778 
BR( 5)= 0.871 
BRC 6)= 0.945 
BRC 7)= 1.010 
BR( 8)= 1.060 
BR( 9)= 1.180 
BR(10)= 1.410 
BR(1i1)= 1.570 
BR(12)= 1.660 
BR(13)= 1.670 
BR(14)= 1.670 
BR(15)= 1.670 
BR(16)= 1.630 
BR(17)= 1.370 
BR(18)= 0.919 
BR(19)= 0.448 
BR(20)= 0.195 
BR(21)= 0.188 
BR(22)= 0.168 
BR(23)= 0.132 
BR(24)= 0.053 
BR(25)= 0.000 


DO 103 M = 1,25 
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XL(M) 
BR(M) 


XL(M)/L 
BR(M)/L 


103 CONTINUE 


101 


PISIM =TSIM/DELT 
ISIM =PISIM 


SIMULATION BEGINS 


DO 100 I=1,ISIM 


CALCULATE THE DRAG FORCE, INTEGRATE THE DRAG OVER THE VEHICLE 


DO 101 K=1,25 
UCF=W-XL(K) *Q 
VEC1(K)=BR(K) *UCF*ABS (UCF ) 
VEC2(K) =BR(K) *UCF *ABS (UCF ) *XL(K) 
CONTINUE 
CALL TRAP(25,VECi ,XL,HEAVE) 
CALL TRAP(25,VEC2,XL,PITCH) 
HEAVE=CDZ*HEAVE 
PITCH=CDZ*PITCH 


COUPLING EQUATIONS 


DWDV=(IY-MQDOT) *(-HEAVE) +(MASS*XG+ZQDOT)*PITCH 
DQDV=(MASS-ZWDOT ) *PITCH+ (MASS*XG+MWDOT) *(-HEAVE) 
CiDV=CIY-MQDOT) *MASS*ZG*Q**2- (MASS *XG+ZQDOT) *MASS*ZG*W*Q 
C2DV=-(MASS-ZWDOT) *MASS*ZG*W*Q+(MASS*XG+MWDOT) *MASS*ZG4Q**2 


DW=DWDV/DV 
DQ=DQDV/DV 
C1=C1DV/DV 
C2=C2DV/DV 


WOOT =A11*W+A12*Q+A13*(XGB*COS (THETA) 
+ZGB*SIN(THETA) )+DW+Ci 

QDOT =A21*W+A22*Q+A23*(XGB*COS (THETA) 
+ZGB*SIN (THETA) )+DQ+C2 

THEDOT=Q 

IF (I.GT.3) GO TO 150 
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150 


INITIAL FIRST ORDER INTEGRATION 


W = W + DELT*WDOT 

Q = Q + DELT*QDOT 
THETA = THETA + DELT*THEDOT 
Fi(I)=WDOT 

F2(1I)=QDOT 


F3(1)=THEDOT 
ADAMS-BASHFORTH INTEGRATION 
F1(4)=A11*W+A12*Q+A13* (XGB*COS (THETA) +ZGB*SIN (THETA) )+DW+C1 


F2(4)=A21*W+A22*0+A23% (XGB*COS (THETA )+ZGB*SIN (THETA) )+DQ+C2 
F3(4)=Q 


W =W +(DELT/24.0)* (55 .0*F1(4)-59 .0*F1(3)+37.0*F1(2) 
-9.0*F1(1)) 

Q =Q +(DELT/24.0)*(55 .0*F2(4)-59 .0*F2(3)+37.0*F2(2) 
-9 .0*F2(1)) 

THETA=THETA+(DELT/24.0)*(55.0*F3(4)-59.0*F3(3)+37 .0*F3(2) 
-9 .0*F3(1)) 


F1(1)=F1(2) 
F1(2)=F1(3) 
F1(3)=F1(4) 
F2(1)=F2(2) 
F2(2)=F2(3) 
F2(3)=F2(4) 
F3(1)=F3(2) 
F3(2)=F3(3) 
F3(3)=F3(4) 


TIME=I*DELT 

J=J5+1 

IF (J.NE.IPRNT) GO TO 100 

ALFA=W 

ALFA=(-ATANCALFA) )*180/PI 

WRITE (20,991) TIME,THETA*180/PI,Q,ALFA 
J=0 
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100 CONTINUE 


991 


STOP 
FORMAT (4E15.5) 
END 


SUBROUTINE TRAP(N,A,B,OUT) 
NUMERICAL INTEGRATION ROUTINE USING THE TRAPEZOIDAL RULE 


DIMENSION A(1),B(1) 

N1i=N-1 

QUT=0.0 

DO 1 I=1,N1 
OUT1=0.5*(ACI)+ACI+1))*(BCI+1)-B(I)) 
QUT =OUT+0UT1 

CONTINUE 

RETURN 

END 
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